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SUMMARY 


The Christensen Theory of a stochastic model for hydrodynamic 
lubrication of rough surfaces is extended to elastohydrodynamic lubrica- 
tion between two rollers. The Grubin-type equation including asperity 
effects in the inlet region is derived. Solutions for the reduced pres- 
sure at the entrance as a function of the ratio of the average nominal 
film thickness to the r.m.s. surface roughness (in terms of standard 
deviation a), have been obtained numerically. Results were obtained for 
purely transverse as well as purely longitudinal surface roughness for 
cases with or without slip. The reduced pressure is shown to decrease 
slightly by considering longitudinal surface roughness. The transverse 
surface roughness, on the other hand, has a slight beneficial effect on 
the average film thickness at the inlet. 

The same approach was used to study the effect of surface roughness 
on lubrication between rigid rollers and lubrication of an infinitely- 
wide slider bearing. Results of these two cases show that the effects 
of surface roughness have the same trend as those found in elastohydro- 
dynamic contacts. 

A comparison is made between the results using the stochastic approach 

and the results using the conventional deterministic method for the inlet 

pressure in a Hertzian contact assuming a sinusoidal roughness. It was 

found that the validity of the stochastic approach depends upon the 

number of wave cycles n within the Hertzian contact. For n larger than 

a critical number, which depends upon the ratio of asperity height to the 

nominal film thickness, 6 /h , the stochastic theory yields the same 

max o 

results as that obtained by the deterministic approach. 
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Using the flow balance concept, the perturbed Reynolds equation, 
which includes a single three-dimensional rigid asperity in one of the 
lubricating surfaces, is derived and solved for the perturbed pressure 
distribution. In addition, the Cheng's numerical scheme , for EHD 
contacts, is modified to Incorporate a single two-dimensional elastic 
asperity, or a waviness profile, on the stationary surface. The perturbed 
pressures obtained by these three different models are compared. Quali- 
tatively, the results for the single 2D elastic asperity and the waviness 
profile model by using Cheng's scheme, are mostly the same. However, 
some results obtained for the single 3D rigid asperity exhibit different 
trends when conq>ared with the single 2D elastic asperity or the waviness 
profile. In the case of the wavlness profile in which the local elastic 
deformation is allowed, the magnitude of the pressure fluctuation ^ , is 
found to increase when the pressure viscosity parameter G, or the ratio 
of the asperity anplitude to the nominal film thickness, 6^^/h^, increases. 
On the other hand, is found to decrease as the magnitude of the Hertzian 
pressure F , or the ratio of the nominal film thickness to the radius of 
the equivalent cylinder, h^/R increases. 
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CHAPTER I 


‘ ^ INTRODUCTION 

In conventional sliding bearings, there usually exists a higji degree 
of conformity between bearing surfaces, and this enables a substantial 
load to be generated by the thin oil film. The performance of these 
bearings can be satisfactorily predicted by solving the Reynolds equation 
for the pressure distribution within the lubricant film^ However, for 
highly loaded concentrated contacts, such as gears, cams and rolling 
contact elements, the lubrication phenomenon cannot be predicted by 
Reynolds equation alone. Local elastic deformation of the solid under 
high pressure becomes influencial in determining the load capacity and 
film thickness of these contacts. The study of lubrication processes 
including the elastic effects is presently known as elastohydrodynamic 
lubrication (EHL). 

To date, the theories for hydrodynamic and elastohydrodynamic lubri- 
cation have reached a very advanced stage. However, these theories are 
mostly based on the assumption that the lubricating surfaces can be 
described by smooth mathematical functions. In reality, surfaces are 
never perfectly smooth in a microscopic scale. In the hydrodynamic 
lubrication regime, the asperity heights of the rough surfaces are much 
smaller than the average lubricant film. Thus, the effect of surface 
roughness on hydrodynamic lubrication, in most cases, can be neglected. 
Hence, the smooth film hydrodynamic lubrication theories provides a very 
satisfactoiy prediction of lubrication performance. In elastohydrodynamic 
lubrication of concentrated contacts , there exist two distinctively dif- 
ferent regimes, the full film and the partial film EHL. In the full film 
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regime, the average nominal film thickness is usually much greater than 
the asperity heights, and, in this regime, the behavior of the contact 
can be predicted quite satisfactorily by smooth-film EHL theories* In 
the partial film regime, the asperity heights are of the same order as 
the average nominal lubricant film. Thus, the effect of surface rough- 
ness in the regime must be considered, 

EHD film thickness has been well accepted as an important bearing 
design parameter. The degree of asperity interactions, and the surface 
distress in the forms of wear, pitting and scuffing are associated with 
the ratio of the average nominal film thickness to the r,m,s, surface 
roughness (in terms of standard deviation a), in EHD contacts. In many 
cases, bearing failures can be attributed to insufficient film thickness 
which leads to asperity contacts. Thus, there is a need to determine 
the surface on the film forming capability in EHD contacts. Therefore, 
the second chapter of this dissertation is focused on the effect of sur- 
face roughness on the average film thickness between lubricated rollers. 
The stochastic theory developed by Christensen [7] is extended to deter- 
mine the surface roughness influence on the inlet film thickness of EHD 
contacts. In addition, the roughness effect on the load capacity in 
rigid rollers and infinitely-wide slider bearing is also studied. 

The third chapter compares the difference between the roughness effect 
and waviness effect on the average film thickness in EHD contacts. It 
also provides some criteria to determine the applicability of the sto- 
chastic theory. 

The pressure profile enables one to predict the stress distribution 
of the lubricated contact. In the fourth and fifth chapters, the effect 
of surface roughness on the pressure distribution are discussed. In 
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Chapter IV, comparisons are made between the effect on the perturbed 
pressure due to a single three-dimensional rigid asperity and due to a 
single two-dimensional elastic asperity within an EHD contact. In 
Chapter V, the surface profile before elastic deformation is assumed to 
be in the form of sinusoidal waviness. The effects of the following 
non-dimensional variables; the Hertzian pressure Pjjjjj the nominal center 
film thickness h^/R, the pressure viscosity parameter G, the number of 
wave cycles within the Hertzian contact, n, and the asperity height h^/R, 
on the magnitude of the pressure fluctuation as well as the perturbed 
pressure profile are studied. 
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CHAPTER II 


, THE EFFECT OF SURFACE ROUGHNESS ON THE- AVERAGE FILM- ■ 

THICKNESS BETWEEN LUBRICATED ROLLERS ■ * 

2,1 INTRODUCTION 

The Inclusion of surface irregularities in lubrication analysis can be traced 
back to [1-3], in which the roughness is modelled as sinusoidal or saw-tooth 

j ■ 

waviness. Subsequently, Tseng and Saibel [63 introduced the stochastic concept 
based on random surface roughness analysis on lubrication. Their method deals 
with surfaces with one dimensional transverse roughness only. The stochastic 
model has been revived by Christensen and his colleagues [7,8,9,16,17] in studying 
the lubrication process between rigid surfaces containing surface roughness 
modelled as ridges oriented transversely or longitudinally. Recently, the effects 
striated roughness on both bearing surfaces have been obtained by Rhow and 
Elrod [18]. 

The effect of surface roughness on film thickness in EHD contact has not 
been fully explored. However, there have been some related work. For instance, 
Fowles [4] studied the EHD lubrication between identical sliding asperities. Lee 
and Cheng [5] have studied the effect of a single asperity on the film and pres- 
sure distribution during its entrance into an elastohydrodynamic contact. The 
load sharing between fluid film and asperity contacts as well as the traction in 
partial EHD contacts have been studied by Tallinn [10], and Thompson and Bocci 
[12], Johnson, Greenwood and Poon [11] have used an approximate analysis to 
ascertain the effect of roughness on the EHD film thickness. They concluded 
that, to a first approximation, the separation between two rough surfaces is very 
close to that calculated by the smooth film theory. Recently, pressure and trac- 
tion rippling inEHD contact of rough surfaces have been calculated by Tallian [19], 
by using Christensen’s stochastic model of hydrodynamic lubrication. 
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In the present analysis, Christensen's approach [7] Is extended to determine 
the surface roughness Influence on the Inlet film thickness of EHD contacts. The 
Grubln-type hydrodynamic equation In the Inlet region for the rough surfaces Is 
derived and solved numerically. Results are compared with smooth-film theories. 
In addition, the load capacity In rigid rollers and In an Inflnltely-wlde bearing 
Is also presented for comparison. 
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2 . 2 GOVERNING E QUATION 


Assuming that the lubricant is isothermal and incompressible and the side- 
leakage is negligible, the one- dimensional Reynolds equation governing the pres- 
sure in an EHD contact is 


Sx \12|j, 9x/ \ ,2 / Bx at 


( 2 . 1 ) 


where is the local total film thickness consisting of the following three parts 


h^ = h+ 6 j ^+62 


( 2 , 2 ) 


In the above, h is the local average film thickness, and 6 ^^, 62 , the roughness 
profile measured from the mean level of surface profiles 1 and 2 (Fig. 2.1). 


2.2.1 Transverse Surface Roughness 

In this case the asperities on both lubricating surfaces are straight ridges 
perpendicular to the direction of rolling. Equation (2.2) becomes 


h^ = h + 6 ^(x - u^t) + 62 (x - U 2 t) 


With the relations. 




36^ 

‘i(“ - - "1 — 


Bt 


(x - U2t) = - ^2 



Eq. (2.1) is simplified to 


B r"? ap "1 “2 

8|x L 12|x ax 2 



(2.3) 


(2.4) 


(2.5) 


( 2 . 6 ) 
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It is assumed here that there are enough numbers of asperities within the Hertzian 

zone such that h can be considered as a constant of time. Let the bracketed term 

in the left-hand side of Eq, (2.6) be denoted by 
3 


m.!!L 2£ 

^ 12yi 3X 2 


u. - u, 

h + 2' " ^*1 ” 


(2.7) 


For Inflnltely-wlde slider bearing and rigid roller bearing, M la expressed as 
above. In elastohydrodynamlc contact, the reduced pressure, q , and viscosity, p,, 
which are respectively 
-dp 


1 - e 


( 2 . 8 ) 






(2.9) 


are Introduced. Then, M In an EHD contact Is 
M 


“T aa.ULlIi 

12p, Bx 2 


'"l " "z 

h + 2 ( 6 j^ - 62 ) 


( 2 . 10 ) 


It is shown in appendix A that M in the case of EHD contact, rigid roller 
bearing, or infinitely-wide slider bearing is a stochastic quantity with a negli- 
gible variance comparing to the variance of the terms on the right hand side. 

Re-arranging Eq, (2.7) and taking expected values on both sides, one obtains 


where 




' { } ■ I { } 


( 2 . 12 ) 


and £(r) Is the probability density distribution of the random variable y. 


since M Is a stochastic quantity with zero (or negligible) variance, M and ^ 

T 

can be considered to be (approximately)stochastlcally Independent quantities. 
Hence 


e 




Then Eq. (2,11) can be re-written as 




¥ 


«i +U2 


h + 




(2.13) 


(2.14) 


where p is now the expected or mean value of p. Substituting Eq. (2,14) into 
Eq, (2,6) and re-arranging, one obtains the stochastic Reynolds equation for 
rigid rollers bearing and inf initely-wlde slider bearing. 



Similarly, the stochastic Reynolds equation for EHD contact can be expressed as 


[JL— la 

dx L' 12 p,^ dx 




(2.16) 


2,2,2 Longitudinal Surface Roughness 

VHien the asperity ridges are parallel to the direction of rolling, 6 ^^ and 62 
are independent upon x, u^, U 2 and t. Again, assuming ^ = 0, then Eq, (2,1) can 
be simplified to 



(2.17) 
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For the cases of rigid rollers and inf initely-wide slider bearings, with longi- 
tudinal surface roughness, it is shown in \_7'] that the stochastic Reynolds 
equation is of the form 


u, + u« 


Similarly, the stochastic Reynolds equation of an EHD contact can be shown to be 


dh 

2 dx 


I? « 

If the roughness distribution is symnetric to zero mean 


(2.19) 


(2,19a) 


2.3 METHOD OP SOLUTION 

2.3,1 Elastohydrodynamic Contacts 

Using Grubin’s approach, it is assumed that the average surface profile, h, 
in the inlet region is governed by the deformation produced by a Hertzian 
elliptical pressure distribution in the contacting region. From [14^, this pro- 
file is given by 

^ ~ ° ^ ~ 0 ] ( 2 . 20 ) 

^O o 

Introducing dimensionless variables Q, X, H, H^, a, U, 62 * ^ ^ defined 

in the Nomenclature, Eq, (2.16), for the transverse roughness becomes. 
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The boundary conditions are: 


dX 


Q * 0 


at X = - 1 


at X = - 


(at H = I ) 


( 2 , 22 ) 

(2.23) 


2.3. l.a Pure. Rolling. Case 

In this case, S = 0, and Eq. (2.21) becomes 


dX 


dQ o 
LdX 


o 1 ^ to 

ryJ dx 


“t 


(2.24) 


Integrating Eq. (2.24) twice with the two boundary conditions, one obtains 

* « H - 1 

« - Olx-l - J, 

o 


(2.25) 


3 ri 1 

where G 2 = H e \“ 3 J 


H 


r 








and 


r . ^ 

i. L» + « (|)] 


* — 

6 = 6 /cr 


(2a26) 


6 * ^2 


— -^—2^2 

a is the composite standard deviation and is defined as o “ ^2 * 

g(6 ) is the roughness height distribution function. 
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Once a and g(6 ) 


are given, Eq. (2; 25) can be evaluated numerically for Q , 


2,3. l.b Rolling and Sliding 

If there is relative sliding between surfaces, S will, not be zero, and the 
last term in the stochastic Reynolds equation, Eq. (2.21), will not necessarily 
vanish. However, if the roughness distribution function of both surface profiles 
are the same, then 



(2.27) 


Using this relation, Eq. (2,21) becomes 


dX LdX 



(2.28) 


which is the same as Eq. (2,24) of the pure rolling case, Q will accordingly be 
the same as that expressed in Eq. (2,25). 

If, on the other hand, one of the contacting surfaces is considered rough 
while the other one smooth, then Eq, (2.21) becomes 



(2.29) 


where the plus sign is for = 0, and minus sign for ^2 ~ 0. 

(See Appendix B for the Fortran IV listing of the numerical analysis.) 
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Defining 


, H 

J. ^ 

4 a 


r 




r” g dt* 

-3 J. h" [i + <|)e*y 


•k ^ k^ 

J..s. _(ft ) 


- <l >^7 


(2.30) 


Eq, (2*29), expressed iti G 2 and G^, becomes 


(2.31) 


Using boundary condition (2.22), one obtains 

- 1 +-22 (' i ') 

G 2 dX 2 G 2 ^ 2 


(2.32) 


Integrating Eq. (2.32) between -09 and -1, one obtains the expression for Q 


“Ml ( T") ^^2*^^ ? J ^ [^4 ■ ^2 (gP 

H 1-00 H -CO H 2 X=-l 


dX 


(2.33) 


Eq, (2.33) can be integrated niimerically for Q for various slide to roll ratios 
from the pure rolling case, S = 0, to the simple sliding case for which u^ = u, 

U 2 = 0, and S = 2, 

For EHD contacts with longitudinal surface roughness, the stochastic Reynolds 
equation following Eq. (2.19) becomes 

2 ^ 




dX 


(2.34) 
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The above equation is valid for any rolling and sliding EHD contacts* Using 

* 

boundary conditions, Eqs. (2,22) and (2,23), the reduced pressure Q at X = -1 
can be integrated as 


* 

Q 




dX 


(2.35) 


2.3,2 Rigid Rollers 

If the elastic deformation of the rollers is neglected, the smooth, average 
surface profile can be approximated by a parabolic profile 



where 


(2.36) 


X = x/R 


ho = the average center film thickness 


Using the same approach as developed in EHD contacts, but with different dimen- 
sionless variables for P, X, X , H, H , a, and o as defined in the 

Nomenclature, the stochastic Reynolds equation for rigid rollers with transverse 
surface roughness can be expressed as 



(2.37) 
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This equation has the same form as Eq, (2.21) of the EHD contacts. However, the 
boundary conditions for the case of rigid rollers will be different from Eqs. (2k22) 
and (2.23). They are given by 


^ - 0 , p - 0 at X - X* 

aJ\ 

P=0 atX*-co 


(2.38) 


Using these boundary conditions, Eq. (2.37) can be readily integrated to yield 




(2,39) 


^ ' H' 


2 X*X 


•/If 

when I is a dummy variable for X, and X and H are determined by imposing 

P(X ) ■ 0. Eq. (2^39) can be Integrated numerically to yield P for the pure rolling 

case (S « 0), simple sliding case (S > 2), as well as rolling and sliding case 

•/f 

(any S). Once P(X) Is found, the dimensionless load W , can be determined by 


* r 
w = - 


("l + "2) 


R 


I 

iJ 


w = f 


P(X)dX 


(2.40) 


For rigid rollers with longitudinal surface roughness, the stochastic 
Reynolds equation is 




(2.41) 


Using boundary conditions, Eq. (2.38), one obtains 

* 


C 7T 


H > H 
1+3 




(2,42) 
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^ ^ 

where X and H are determined by the condition P(X ) = 0, The dimensionless 
pressure and load can be determined in the same manner as that for rigid rollers 
with transverse surface roughness. 


2,3.3 Inf ini tely "Wide Slider 

For an infinitely-wide slider, the smooth, average^ surface profile^ can be 
represented by 


H(X) = = 1 + (1 - X) (2.43) 

o o 

where 

h = h . = minimum film thickness at the exit of the slider 

o min 

m = slope of the slider 
Z = length of the slider 


Introducing dimensionless variables P, X, H and 6 as defined in the Nomenclature, 

one obtains the stochastic Reynolds equation for an infinitely-wide slider with 
transverse surface roughness, 



The boundary conditions for Eq, (2.44) are 
P(0) = P(l) = 0 

For both surfaces having the same roughness characteristics. 



(2.45) 


(2.46) 
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Eq. (2.44) becomes 

^ riE ■ is 

dX LdX G 2 J “ dX 


Integrating twice, one obtains 



(2.47) 


(2.48) 


where C Is evaluated by the boundary condition P(l) « 0, For one surface rough, 
and the opposing surface smooth, £q« (2,44) takes the form 


dX LdX 



dX 



(2.49) 


where the plus sign represents the case of a smooth surface sliding against a 
stationary rough surface, and the negative sign implies the rough sliding against 
smooth surface^ Eq. (2.49), expressed in terms of G 2 and 6^, becomes 


is 4. “ i- 

dX LdX G2J “ dX * dX ^2/ 


(2.50) 


Integration of the above equation yields 


r G, r* G. pX G, 

^d|±aj ^d? - cf ^d5 

o H o H o H 


(2.51) 


o H 

where C is determined by the boundary condition, P(l) * 0. 

For the inf initely- wide slider with longitudinal surface roughness, the 
stochastic Reynolds equation is 



(2.52) 
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which, after integrating twice, yields 



(2.53) 


with C determined by P(l) 

•' -i 

by 


r I - 
* p-*- 

W = J P(X)dX 
o 


0, In all the above cases, the load can be evaluated 


(2.54) 


2,3.4 Roughness Distribution Function 

The roughness distribution function employed in this paper is the same as 
that used by Christensen namely, - 


g(6*) = { 


35 

96 9 ^ 


•k * * 

- 6 < 6 < 6 

max max. 


Elsewhere 


(2.55) 


6 =3 

max. 


(2.56) 


This polynomial distribution function is an approximation to Gaussian distribution. 
The reason using this polynomial function is that the roughness height distribu- 
tion function of many engineering Surfaces is very close to Gaussian C20], and 
that a Gaussian distribution always implies a finite probability of having asperities 
of very large sizes which are very unlikely in practice. 
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2.4 DISCUSSION OF RESULTS 


2.4.1 EHD Contacts 

The effect of surface roughness on the pressure generation at the inlet of 

EHD contacts can be presented conveniently by using a quantity K^, defined as the 

He 

ratio of the inlet pressure calculated from the stochastic theory, Q , to that 

K. 

* 

calculated from the Grub in s smooth-film theory, Q . Thus, 

Kj - Q*/^ (2.57) 

In Fig. (2.4), the ratio is plotted against a surface roughness parameter, 

6 , for the following four cases 

max ® 

1) pure rolling with transverse surface roughness 

2) pure rolling with longitudinal surface roughness 

3) rolling and sliding with S = 0,2 and 6^ s= 0 (smooth surface is faster) 

4) rolling and sliding with S = 0,2 and ^2 ” ^ (rough surface is faster) 

The dimensionless load and film thickness for the above cases are W = 3 x 10*^ and 



These curves are obtained by changing the magnitude of 6 .It is seen that, 

max 

for pure rolling with longitudinal surface roughness, is reduced very slightly 

due to surface roughness effects. Even for 6 as high as 0,99, the 

max 

reduction is only about 7.5%, 

Contrast to the effect of longitudinal roughness, Izhe transverse roughness 
has a much more pronounced effect on the integrated pressure. It tends to increase 
the dimensionless reduced pressure and hence also tends to increase the load 
capacity as increases. For pure rolling, the transverse roughness causes 

an increase in Q from 7% to 30% as 6 is increased from 0.6 to 0.99. The 

effect of relative sliding between a smooth surface and a rough surface in EHD 
contacts is shown in the two curves for S = 0,2, Even for such a small slip, there 
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is already noticeable departure from the pure rolling case. For the case \diere 

the smooth surface. is faster, there is an additional pumping effect, compared to 

the pure rolling case. The reverse is true if the roughness surface is faster. 

The effect of is studied in Fig. (2,6). It is shown that the three curves 

for 5 X lO"^, 9 x lO"^ almost coincide with one another. This indicates 

that the roughness effect on EHL is almost entirely independent upon H • 

Fig. (2.6) shows the integrated value of Q against for the condition of pure 

rolling with transverse roughness, for different ratio of 6 ranging from 

max 

0.0 (smooth film theory) to 0.99. It is interesting to note that the curves are 

4c 

parallel straight lines. Again it is readily seen that the magnitude of Q depends 

on the ratio of 6 ^ For 6 =0,99, there is approximately a 20% gain in 

max Diax 

the mean film thickness over that based on smooth film theory. For 6 = 0.9, 

max 

and 0.6, the gain in mean film thickness is about 15% and 57o respectively. 

In the case of single sliding of an EHD contact, i.e. S = 2, elastic deforma- 
tion of asperities begins to be significant at X = -1. At this position, the 

4c 

values of o, the r.m.s. roughness amplitude and g(6 ) the asperity distribution 
function, will no longer be the same as those of the undeformed asperities. There- 
fore Eq, (2,33) cannot be applied under this situation. One should notice that 
Eq, (2.33) is valid for rolling and sliding case, only when the elastic deformation 
of asperities near X * -^1 can be assumed to be negligibly small. This occurs 
only when S is small, 

2.4.2 Rigid Rollers 

The effect of surface roughness on the dimensionless load of rigid rollers 
is presented in terms of which is a quantity defined as the ratio of the 

4c 

dimensionless load from the stochastic theory, to that of the smooth film 
4c 

theory, Wg. Thus 


21 


(2.58) 


r * , * 

Vs 

Results of K_ versus 6 for the following six cases are shown in Fig. (2L.7): 

K max 

1. sinq)le sliding of a smooth roller against a stationary rough roller 
with transverse surface roughness, S = 2, 

2. rolling and sliding with transverse roughness; S = 0.2, 6^ ^ 0 (smooth 
surface is faster), 

3. pure rolling with transverse roughness, 

4. rolling and sliding with transverse roughness; S = 0.2, 6^ - 0 (rough 
surface is faster), 

5. simple sliding of a rough roller against a stationary smooth roller 
with transverse roughness; S = 2, 

6. longitudinal surface roughness. 

-4 

The center film thickness for the above cases is h^/R = 5 x 10 . 

o 

Similar to elastic rollers, the effect of surface roughness on rigid rollers 

with longitudinal roughness is quite small. For *3 ^ 0.99, the reduction of 

max K. 

is only about 57.. 

For pure rolling with transverse roughness, is increased by about 167o 

when "3 *=0.99. For rolling and sliding with S = 0.2 and smooth surface 

max 

moving faster, there is an increase in load carrying capacity, compared to the 

pure rolling case. The reverse is true when the rough surface is faster in the 

rolling and sliding case. When a smooth roller is sliding against a stationary 

rough roller, there is a substantial increase in K_. When "6 = 0.6 and 0.99, 

max 

the gain in are 8.5% and 39% respectively. However, when a rough roller is 

sliding against a stationary smooth roller, is almost unaffected. When 6 

K max 

is greater than 0.66, begins to decrease slightly after a steady Increase. 
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This small dip in is suspected to be caused by using the above mentioned poly- 

nominal function as the surface roughness distribution function. When a sinusoidal 

distribution function [17] which is defined as, 

1 


«(6 5 » { 


( 2 . 59 ) 


/ 9-6 “ - 3 < 6 <3 

0 Elsewhere 

is employed, the results which are not plotted here show that the dip disappears 

and that K_ increases with the increase of "5 

K max 


2,4.3 'Infinitely-Wide Slider Bearing 

In Fig. (2.8), W is plotted against ^ov the following four cases: 

1. simple sliding of a smooth surface against a rough surface with trans- 
verse roughness, 

2. both surfaces having the same roughness distribution functions with 
transverse roughness, 

3. simple sliding of a rough surface against a smooth one, with transverse 
roughness , 

4. longitudinal surface roughness. 

Qualitatively, the roughness effect on a slider checks very well with that on 
rollers. When the roughness direction is longitudinally oriented, the load carry- 
ing capacity of the oil film is reduced.. 

In the case of transverse roughness, the effect on load capacity is always 
beneficial. For both surfaces having the same kind of roughness distribution 

function, W is increased by 667, approximately, when 6 is 0.99. For a smooth 

max 

surface sliding against a rough one, the gain in W is even larger. W is 
increased by about 1297o when 6^^^ is 0.99. For a rough surface sliding against 
a smooth one,W Is only increased slightly. 
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The authors’ results agree very well to those obtained by Rhow and Elrod ClS], 

^o have studied the effects of two-sided straited roughness on the load-carrying 

capacity of an infinitely wide slider bearing. The only exception is that the 

authors’ results show a small dip in W .when 6 is larger than 0,9. This 

max 

small dip in W is caused by employing the polynomial function as the surface 
roughness distribution. When the sinusoidal distribution is used, the dip in W 
disappears. 
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2.5 


CONCLUSIONS 


1, Based on Christensen's stochastic model of hydrodynamic lubrication, a 
Grubln type elastohydrodynamic analysis at the inlet of a Hertzian contact 
indicates that surface roughness can have a noticeable effect on the level 
of mean film thickness between EHD contacts. 

2, For longitudinal surface roughness, ridges parallel to the direction of 
rolling, the present analysis predicts an inlet pressure or inlet film 
thickness smaller than that predicted by the smooth-film EHD theory. 

3. For transverse surface roughness, ridges perpendicular to the direction of 

rolling, the inlet mean film thickness level is increased noticeably due 

to the additional pumping by transverse ridges. The level of increase is 

mainly a function of ratio of the maximum ridge height to the 

mean film thickness at the inlet and is not sensitive to other operating 

parameters. For 6 /h approaching unity, which corresponds to h /a * 3 
max o o 

for 6 = 3a, one can expect an increase of 257o in mean film thickness 

max * ^ 

compared to the smooth-film EHD film thickness for pure rolling. Results 
for small slide to roll ratios indicates that the case of smooth sliding 
over rough gives further enhancement in mean film thickness, whereas the 
case of rough sliding over smooth yields a slight reduction in mean film 
thickness comparing to pure rolling. For high slide to roll ratios, it 
was found that local elastohydrodynamic effect due to local pressure fluc- 
tuations will become significant, and the present analysis based on the 
Gnibin approach will become invalid. 

4. For rigid rollers and inf initely-wide slider bearings, load capacities cal- 
culated for rough surfaces show trends similar to those found in EHD contacts. 
However, the effects for inf initely-wide slider bearings are much stronger 
than that for rigid rollers. 
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TRANSVERSE ROUGHNESS 


U2 



LONGITUDINAL ROUGHNESS 


Figure 2-1 EHD Contacts Between Two Rough Surfaces 
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TRANSVERSE ROUGHNESS 

Figure 2-2 Rigid Rollers With Rough Surfaces 
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TRANSVERSE ROUGHNESS 


Figure 2-3 An Inf inltely-Wide Slider Bearing with Rough Surfaces 
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E H D CONTACTS 
POLY. D I ST. 

W = 3 X 10® Hq= 10 ° 



Fig, 2-4 The Effect of Transverse and Longitudinal Roughness on the 

Ratio of the Reduced Pjes^ure for Rough Surfaces to that for 

Smooth Surfaces, K_= Contact 

E K 'D 
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RIGID ROLLERS 
POLY. DIST. 
hp/R - 5 X 10'* 



* 

Fig. 2-7 The Variation of the Normalized Load Ratio, K = ^ /W^, with 
the Roughness Height, § ^ /h , for Rigid Rollers 
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0.99 


^max 


Fig, 2-8 The Variation of the Nor alized Load, W = — ^ ^ -^-w, 

with the Roughness Height, for an Inf initely-wide 

Slider Bearing 
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CHAPTER III 


WAVINESS AND ROUGHNESS IN EIJISTOHYDRODYNAMIC LUBRICATION 

3.1 INTRODUCTION 

A stochastic theory for elas tohydrodynamic lubrication of contact with two- 
sided roughness has been developed in Chapter II. The basic requirement of this 
theory is that the roughness pattern must be very dense within the contact zone. 

In other words, the largest wavelength in the roughness spectrum must be small 
compared to the contact width. If the largest wavelength is of the same order of 
the contact width, surface roughness becomes surface wavit\ess. The stochastic 
theory may be invalid in this region. 

In order to ascertain the conditions under which the stochastic theory becomes 
valid the reduced pressure based on the deterministic approach using a sinusoidal 
profile and that calculate^ from the stochastic theory using a surface roughness 
distribution equivalent to the sinusoidal profile is compared. The comparison 
is only made for the transverse roughness since the effect of longitudinal rough- 
ness is usually negligibly small comparing to the effect of transverse roughness. 

3 . 2 GOVE RNING E QUAT ION 

In comparing the effect of waviness and roughness on the reduced pressure, 
the surface profile is assumed to be in the form of sinusoidal waviness. The 
reduced pressure is then solved by the deterministic approach for the waviness 
case. At the same time, a density distribution function equivalent to the chosen 
sinusoidal wave is evaluated. Using this equivalent roughness distribution func- 
tion, the reduced pressure for the roughness case is then solved by the stochastic 
approach as illustrated in Chapter II. Then, the effect of surface waviness and 
roughness will be compared for the pure rolling case. 
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3,2,1 The Waviness Case 

In the waviness case, the one-dimensional Reynolds equation governing the 


{pressure in an EHD contact of an isothermal and incompressible lubricant is 
Bx \12a 3xy 2 Bx ^Bt 


where h^ is the local total film thickness consisting of the following three parts 
h, 6^, #2 

where h = the nominal smooth part of the average film thickness 

^1*^2 ” roughness amplitude measured from the mean level of surfaces 


1 and 2 

H [(2n^TT)(x - u^t)] 

^2 “ ®max2 


(3.2) 

(3.3) 


Uit 


= U^t 


It is assumed that the asperities on both surfaces are straight ridges perpendicular 
to the direction of rolling. With the relations 


B6^ 

at” “ ■ '^1 ajT 


(3.4) 


362 562 


(3.5) 


Eq, (3,1) is simplified as 



(3.6) 


It is further assumed that there are enough number of asperities within the 
Hertzian contact zone such that h can be considered as a constant of time. 
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i.e. 




(3.7) 


Hence, in the case of pure rolling, one obtains 

.3 


dx \ 12 g, dx/ " dx 


(3.8) 


with dimensionless variables Q^, X, H, H^, a, U, 6 j^, 62 * ^ defined in the 

Nomenclature, Eq. (3.8) is transformed to 


_ 1 dH 

dX dX J 


T dX 
H 

o 

where H, the average surface profile in the inlet region is governed by 
h - h 


H - 1 


o o 

The boundary conditions are 


= ^[| x 1 a / x 2 - 1 - in{\x\ +Jx^- 1 ) ] 


dX 


at X = - 1 


(3.9) 


(3.10) 


(3.11) 


Q„“0 


at X = - « 


(3.12) 


Integrating Eq. (3.10) twice with these two boundary conditions, one obtains 

rH - 1 . 

(3.13) 

«; ■ 

In the pure rolling case 


^ “ Owl = J (rr-r) 

^ ^lx=-l -® H H„ ' 


9j^=e2=9=«t 


In addition, it is assumed that 


(3.14) 


n^ ^2 ^ 


max 1 max 2 max 


(3.15) 

(3.16) 
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For given W, n and 0, Eq, (3,13) is solved numerically by Simpson's 

integration routine. 


3.2,2 The Roughness Case 

The probability density function corresponding to a sinusoidal wave is 1^173 

1 


TT 9 


- 3 < 6 <3 


f(6 ) 


(3.17) 


elsewhere 


From Chapter II, the corresponding reduced pressure at the inlet is 
-1 

.dX 


* r ^ \ 

On “ J \““2 3 ^ / ^2 

o 


(3.18) 


where 


3 


d6 


(3.19) 


* * 


0^ = Q calculated by stochastic theory. 

★ 

Hence Q., is evaluated for different values of 6 

^ max 


3.3 DISCUSSION OF RESULTS 

Though the following examples are only connected with the pressure generation 
at the inlet of elastohydrodynamic contacts, yet, qualitatively, the general trends 
of the results will be relevant to other t 3 q>es of bearing and surface irregularity. 
The comparison between waviness and roughness can be presented conveniently 

by using a quantity which is defined as 

* * 

0^ “ Qd 

V « ^ ^ (3.20) 
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where and stand for the inlet reduced pressure calculated from the waviness 

model and the roughness model respectively. Hence, by definition, is the 

fractional deviation of the inlet reduced pressure of the Waviness model from 

that of the corresponding roughness model. ^ The results of such a comparison are 

shown in Fig. 3.1 to Fig. 3.3. The dimensionless load and film thickness for 

these examples are W = 3 x lO"^ and H = h /R = 10 while the ratio 6 /h in 
^ o o * max o 

these three figures are 0.3, 0.45 and 0.6 respectively. The phase angles chosen 
are 0, tt/2, tt and - rr/2 while the n- values are integers. 

It is readily seen that heavily depends upon the phase angle 0 for 
small n. Particularly, for 0 = rr/2 and - tt/2, the magnitude of even changes 
signs. However, the effect of phase angle on decreases rapidly as n increases. 
Furthermore, for larger 6^^^ which means more pronounced asperity interaction, Kj^ 
is larger for the same n. For smaller ^ is smaller. The approximate values 

of at n = 5 and n = 10 for the extreme cases of 0 = - rr/2 and rr/2 are listed as 
follows 



n = 5 

n = 10 

^max^^o 

0.3 

0.45 

0.6 

0.3 

0.45 

0.6 


±3.57.1 

± 67. 

±107. 

±0.57. 

± 17. 

±1.57. 


These results provide a better understanding to the statistical roughness 
theory. First, they show that the effect of phase angle vanishes with increasing 
n. Second, the effect of n becomes less important when n is greater than some 
critical value for a given discrepancy between the waviness model and the rough- 

ness-model becomes poorer as 6 /h increases. It is quite evident that n and 

^ max o ^ 

6 /h are both important parameters that determine the validity of the statistical 
max o 

theory for roughness surfaces. For the particular numerical example used in these 
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II 


calculations. It Is found that for n > 10, waviness is equivalent to roughness and 
that the stochastic theory holds. Even for n * 5, one can still apply the sto- 
chastic theory with reasonable accuracy. 
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3.4 CONCLUSIONS 


For a given W and 

1. 0, n and 6 /h are the parameters to determine the deviation of 
the stochastic, roughness made from the deterministic waviness 
model . 

2. The effect of phase angle diminishes with increasing n. 

3. The effect of n vanishes as n becomes large, and the stochastic 

theory for roughness surface is proved to be valid as n approaches 
a critical value depending on for a given W and . 
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CHAPTER IV 


PRESSURE PERTURBATION IN EHD CONTACTS 
DUE TO AN ELLIPSOIDAL ASPERITY 


4.1 INTRODUCTION 

Recently, there has been a growing interest in the effect of surface 
roughness on the bearing performance in thin film lubrication. In Chapter 
II, the stochastic theory is used to study the effect of surface roughness 
on the average EHD film thickness and the integrated pressure at the inlet 
of the lubricated Hertzian contacts. However, the stochastic theory is in- 
capable of predicting any detailed local perturbations in pressure or de- 
formation caused by the asperities. It was recently pointed out by Tallinn 
[19] that the pressure ripples can rise to a very high level in rolling and 
sliding EHD line contacts. These ripples are very likely one of the chief 
attributing factors to contact fatigue. 

In the last few years, there has been considerable interest in the 
basic event involving a single asperity entering an EHD contact [S] or the 
encounter betoeen two identical asperities ^4]. The work in this chapter 
is aimed towards gaining further understanding of the effect of a single 
asperity on pressure distribution in a line EHD contact. Special attention 
is given to the three-dimensional aspect of the asperity which is asstuned 
to be ellipsoidal at the tip. The effects of ellipticity (aspect ratio) on 
the double amplitude of pressure fluctuations under various rolling and 
sliding condition is examined in detail. 
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4o2 MATHEMATICAL ANALYSIS 


The present analysis consists of two parts. The first part studies the pres- 
sure fluctuations due to a single three-dimensional ellipsoidal asperity at the 
inlet region of an EHD contact assuming that the asperity shape is unaffected by 
the perturbed pressures. These pressure fluctuations are determined by solving a 
perturbed Reynolds equation, in which the unperturbed pressure profile is obtained 
by using a line contact EHD analysis [22], Results are presented as the perturbed 

pressure profile, as well as the amplitude of the pressure ripple, A s as a 

s 

function of ellipticity ratio, y, maximum Hertzian pressure, nominal EHD film 

thickness h^/R, asperity size, b, asperity height, c^/h^, pressure viscosity coef- 
ficient, G, slide to roll ratio, S, and the position of the asperity center X^, 

In the second part, the line contact EHD analysis [22] is modified to include 
a two dimensional asperity ridge on the stationary side of the lubricated contacts. 
In this approach, the elastic deformation of the asperity is included. Results 
which are presented as the double amplitude of the perturbed pressure as a function of 
^l^^o’ compared with those obtained for the three 

dimensional ellipsoidal asperity with large ellipticity ratio for the case of 
simple sliding between a smooth surface and a stationary asperity (S=2)o 


4,2,1 Geometrical Configuration 

The contact between two cylinders as sho\m in Fig. (4,1a) can be described by 

an equivalent cylinder near a flat surface as sho^-m in Fig. (4.1b), As the contact 

width is very small compared to the dimension of the cylinder, the film thickness 

for a rigid cylinder, h , without the asperity is 

§ 


h 

g 


h 

o 



(4ol) 


where x = coordinate along the film 
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R^R^ = radius of rollers 1 and 2, 

h = undeformed center film thickness 
o 

With elastic deformation, the film thickness profile, still without the asperity 
[223 becomes 

2 ^f 

»1 “ J. toi^p^«)d5 (4.2) 

where h^ = smooth-film thickness 

ho = center film thickness at x=0 



= Young's modulus for rollers 1 and 2 
= Poisson's ratio of rollers 1 and 2 
^ = dummy variable for x 

P-j^(5) = smooth- film pressure profile 

Referring to Fig. (4.2 ) the height of a three-dimensional ellipsoidal asperity 
can be x^ritten as 

6 = 6^cos (4o3) 

As the contact width is very small compared to the radius of the cylinder, T\ is 
very small and 

cos T] 1 (4*^) 

Thus the asperity height function of a three-dimensional asperity can be approxi- 
mately \7Titten as 
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( 4 . 5 ) 


,x - x„. 2 




6(x,y) “ 6j^ 


'0 elsewhere 

Similarly, the asperity height function of a two dimensional asperity ridge can 
be expressed as 


I (x) =- I 


X - x^ 


f4(”b^) " ^ 1 


I 

elsewhere 


( 4 . 6 ) 


The total local film thickness, h^, at any point under the surface 
asperity is 




(4.7) 


4.2.2 Governing Equations 
4. 2.2.1 The Smooth-Film Case 

Referring to [22^, the two coupled equations governing the pressure and film 
distributions in an elastohydrodynamic line contact between two rollers with iso^ 
thermal and incompressible lubricant are; 

The Reynolds equation 

dpi ,h - h . 

— = 6u(u^+ U2) (-3 ) ( 4 .; 

where p^, h^, h and x are already defined 
p, = viscosity 

^l’^2 ^ velocity of rollers 1 and 2 


and the film thickness profile as described in Eq. (4.2). 
Eqs, (4.8) and (4.2) become 



In non-dimensional form. 


( 4 . 9 ) 
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(4.10) 



I » |/b, = x^/b. 

Using the same numerical scheme as developed in [22], these two equations are solved 
by the Newton-Ralphson method* P^(X), H^(X) and U are solved for a given set of 
h /R, P , and G* These smooth-film results are used as the inputs to the perturbed 
Reynolds equation described later. 


4.2*2*2 The Discretized Reynolds Equation 

It is assumed here that the viscosity is an exponential function of the pres- 
sure with a pressure viscosity coefficient or, i.e. 


* M's 




(4.11) 


In numerical form the steady Reynolds equation for an incoirpressible lubricant 
can be expressed by dividing the contact zone into small grids with irregular 
spacings* Referring to Fig. (4.3) and applying the principle of conservation of 
mass, one obtains for the ith grid 


where 


®l" ** ^restored 

“l “ “ \12jr ,,, 2 ^P^l^i-1/2 


(4.12) 


V '^2 


i-1/2 
-aPi 3p 




/ 3 -orPi op^N 


U,+ u„ 


i-1/2 


2 ^*^1^ i-1/2 


(4.13) 
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“2 = 



■QfP]^ 


^1+1/2 ° 2 <'>l>l+l/2 


m 


restored bt l 


= p<^>ii 2 xj = ° 


Combining Eqs. (4.12) to (4.15) and re-arranging, one obtains 


1 r(h3 /""i - (' p'. ] 

"" 2 [^’^l^i+1/2 ■ ^’^l^i-l/2j 


(4.14) 


(4.15) 


(4.16) 


Thus p^, the smooth- film pressure profile and h^, the smooth- film thickness are 
functions of x only. 


4. 2. 2, 3 The Perturbed Reynolds Equation 

The pressure distribution can be considered as the sum of the smooth-film 
pressure, p^ and the perturbed pressure, 0. Define 

p = p^+ 0 (4.17) 


Since p^ is a function of x only, the derivatives of p are: 

+M 

Bx Bx Bx 


(4.18) 


|£ = M (4.19) 

9y ay 

For a two-dimensional flow field, the principle of conservation of mass as applied 
to the (i,j)th grid as shown in Fig. (4.4) yields 
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(4.20) 


[ ^2 ”2> + [ ^ (“3- “4> 

= 2^{(ph) 

at l'-PVi-l/ 2 ,j L 2 J 2 

4. ri, ^ r <^y)i-i-*- <^y>n 

+ (pVi+l/2,j L ^“2 2 J 


where m 


pK 


, U ,+ Uo 


r "'“T ap . / I “2^ , I 

1 = T5p: ^ + p l-”2 — ) V, 




, ph^ .o,(p^+|) u^+ U2. , 

= V 12T ® + p — ) (^1+ «>;. ,,, . 

s 1 - 1 / 2 , J 

It is assumed that the value of of0 is much smaller than unity such that e 
be linearized as 

= (1 - 

by using Taylor's series expansion. Neglecting second order term, Eq. (4.21) 
be re-written as 


J- Ph^ -orPi 30 oPlx /-i- -2^ . 


9p, 


,u,+ u„ 




Likewise, m 2 , m^, and m^ can be shown to be 


“2 






ph -op, zap. 


,u,+ u„ 


1 ^ 4 . ri 2> . .,'1 


I+1/2.J 


. . /“"i 24 

LIV. »Ji.3-l/2 


m, = 




Ll2y 




i» j+1/2 


(4.21) 
can 

(4.22) 
can 

(4.23) 

(4.24) 

(4.25) 

(4.26) 
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Substituting Eqs, (4.23) to (4.26) into the left-hand side of Eq. (4.20), one 
obtains 








3 ^ 3 \ "“Pi ^Pl> 


( Vli) } 


, p (it-) [(-T 1^) 


i, j+1/2 


' ^ ^y'i,j-l/2-' 


(4.27) 


With the relations 

f^h, 


(4.28) 


and 


^ 6(x - u„t) = - u„ 


St 


2 9x 


(4.29) 


the right-hand side of Eq. (4.20) can be re-written as 


(4.30) 


Equating Eqs, (4.27) and (4.30) and simplifying, one obtains 
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(.y)j (.y) J { r (h? - (hj II) ] 

j-1 jj L ^ \ T ‘*'^^i+l/2,j ^ ^ ‘*’"^i-l/2,j-‘ 


[-/, 3 "“Pi ^Pl-^ ' / 3 "“Pi *^Pn 1 \ 

LVS ® ■" V\ ® dTj. , J J 

1+1/2, j i-l/2,j 






i+l/2,j 


r« - >.?) ^]. ^ ^.(v u,) 

1 - 1 / 2 , J 


Using the central difference approximation for the pressure gradients, Eq« 
is discretized to 

L(Ay)j_i+ (Ay)jJ { 


j 


±il ' / ’^i+l,j~ '^i.j 
J V (AX). 




3 "“Pi \ / 3 "“Pin 

K ^ ) + (h;; e ^ ^ 

2 J 




/ 3 "“Pi N / 3 "“Pi N 

(h e <>«) +^>**Th '^). ., r<Pl>. - 

. r ^+i »3 iij. ! r 1+1 1 1 

L 2 J L (Ax). J 

(hj e ^0,0) +(h^e ^ a^) (Pi)/ (Pi> . 

-[ ^^^-2 [■' ) 




(4.31) 

(4.31) 


(4.32) 
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(h! + (hi , \ 

^ <h^) J { [1^^ 




("? 


i) + (hi .■"!') (^ , \ 

Ua ^ ^ ^. i - i - | W.r ^ 1 . 1 - 1/ 1 

2 J to>].i J 


[(Ay)j_i+ (Ay)j] { 

[(»? - -?) e'"^^ 

+ 

[(-^1 - »S 

Li 

[ (Pi> - (Pi> ] 

^ i+l i"" 


2 



(Ax). 


[(h^ - hj) ® ^] + (^1 " *^ t ) ® T 

^ -‘i.1 ^ ^ -^1-1. i r ^1 ^ 1-11 

2 L J 


, up (!i±i 4 ^) 


(4.32) 


After re-arranglng and the following non-dimensional variables are introduced: 


^ = Ax/b', IX “ ^y/d', Y = a'/b', Hj,“ 


x/b. 


X, = x-/b, 




Y = y/a' , b = b'/b, 

2.2 

, n / 

E' 


U„ 


M<S^V "2> 
2E'R 


hj^/h^, Cj^/Hq , 


Pr Pi^Prz* ^ “ ^''Prz* 




O' = aPjj2. 


,x - X-x 2 


“o-'’o'h. Sro" *« V®o -h. 

X - X,x 2 


a - ./h . cp.X(^) . (i.) - i] - c, [(V^) . (x>^. i] 


the governing equation becomes: 

[ ^V ^ l { (4 (4 , , 


+ a (h^ e [ (p ) - (p ) ] ]■ ? ^ 

^ i-l,r i i-1’^ 


(4.33) 
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1 r (&X)^_3^+ 3 

*'■ 2 1 . (AY) . 1 JLv®T 


),,/(«? - 

.•1 3-»J . ! ±ZJl2 


+ {-[tov)j.i+teY).][ 


,(h? («? 




. - ^„3 '°^l\ 

+ v“t ® ; . . 

3 ->J 


- ij[(iX>i.i+ (i«)iJ ■ 


<p.>.„-<pi), ,, <n),-<pi>„ 


- /..3 “t'lN " i x-11 

- “ v«T ^ y . . — (AX). - J 

1 >J 


(■4 


/ 3 ""Pin 

' +K ® ). . 

i.i+i i.j 

(AY)/ 


i , 1 k»J 






1 p (A^) • 1 + (AX).-| p/ o ”(VP]^\ / 3 "] 




_ / o ”^Pi\ r '1 ^ 

- " (®T ® i ^i+i,j 


i+l, j 


i+1 i 


. <m"\ ( { [(hJ - 4) . 


- H?) [ 


^Pp..,' ^Pp. 


^] - {K - H?) 


(AX)^ 


(Pi) “ (Pi) 


[(-J - »?) ^ , } [^ferr^J ^ ^ 

1 - 1 , J 1-1 




(4.33) 
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4,2,3 Method of Solution 


In Eq. (4*33), y ^i+l,j unknowns. 

and considered as known are determined by the method outlined in [22^ for the 
smooth lubricated contacts. Thus, Eq. (4.33) can be re-written in the following 
matrix form: 


[*i] * ['ll {’j + ['ll { »i+i} - 

where [A^] and [c^] are M x M diagonal square matrix- [B^] is a M x M tri-diagonal 
square matrix- ^^i+1^ M x 1 column matrix^M and N are 

the number of grids used in the x and y directions. 

In prescribing the boundary conditions for Eq, (4.34), it is necessary to 
assume that the effect of the asperity on the pressure distribution at |X-x^/b*| S: 5 
or |y/a*| ^ 5 is negligible. This assumption justifies the prescription of the 
contact. Thus the boundary conditions are $ = 0 along X-x^/b* = ± 5 and along 
y/a’ = ± 5, or 


{»i} ■ w * " 

(4.35) 

f . . = 0 for j = M 

^ > J 


Along 3^=0 or j=l, the flow is considered to be symmetric, and this gives the follow- 
ing additional boundary condition. 




for j=l 


(4.36) 


Eq. (4,34) is solved numerically by the Columnwise Matrix Inversion Method [25]. 


4.3 DEFORMED ASPERITY ON THE STATIONARY SURFACE 

In the case of simple sliding of a smooth surface against a stationary asperity, 
the non-dimensional Reynolds equation and elasticity equation are respectively 
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(4.37) 


dP _ 48^ „ - 

dX '' 2 J ^ \ 3 / 

“t 

H = H + — (4.38) 

^ ^ h 

o 

In this case, 6 is a function of x only. Thus, ^6/bt = 0. Therefore these coupled 
equations can be solved siinultaneously using the same numerical scheme as described 
in [22']. ^The results are compared with those obtained by solving the perturbed 
Reynolds equation. 


4.4 DISCUSSIONS OF RESULTS 


The results of the perturbed pressure $ are presented as , which is defined 
as 


s max min/ 


at j=l 


(4.39) 


where the position j«l is at the centerline of the asperity along the sliding 

direction. At this position, the effect of asperity on $ is most severe. Since 

$ is a function of P , G, h /R, b, Yj c- /h , x„, and the slide to roll ratios, 

Hz o i O J 


the effects of these variables on ^ are studied separately. 

s 


(1) Effect of c/R 

With P = 0.003, G = 100, X« = -0.5, b = 1/32, c /h = 0.3, the results of 

HZ j 1 o 

h^/R versus are shown in Fig. (4.6). In the case of 2-D elastic asperity, the 

— 6 —6 

magnitude of A increases as h /R increases from 1.0 x 10 to about h /R = 5 x 10 
so o 

(approximate) after which the magnitude of decreases with further increase in 
h^/R. This phenomenon is attributable to the local elastic effect of the asperity. 

In general, one might imagine that the effect of roughness on f should be much more 
severe as the ratio h^/R becomes smaller. However, in an elastohydrodynamic contact. 
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the elastic effect becomes more significant as h becomes smaller. Thus, for a 

fixed ratio of c^/h^, the elastic effect tends to flatten out the asperity and 

produces a trend which shows a decrease in as h /R decreases. On the other hand, 

s 

for the 3D rigid asperity analysis the elastic effect is not accounted; therefore, 

the magnitude of ^ consistently increases with a decrease in h /R, 

8 o 

Based on the results shown in Fig. 4.6, it appears that further comparisons 
between the 2D-elastic asperity and 3D-rigid asperity results are more meaningful 
for values of h^/R greater than 1,0 x 10*^, since below this value no good agree- 
ment is expected. 

(2) Effect of P 

, rl2 

In hydrodjmamic lubrication, with a rigid asperity, it is expected that the 

effect of asperity on will become more severe as the load increases. However, for 

s 

EHD contacts, this simple trend is only expected to be true if the load parameter 
P„ is sufficiently small. For very heavy loads, the elastic effect will iron out 

HZ 

the asperity and this may decrease the magnitude of ^ . Such trend is demonstrated 

s 

in Fig. (4.7), In the 2D-elastic asperity curve, the magnitude of A_ increases 

s 

with P„ until P « 0.0045. Beyond this point, the magnitude of flattens out 
HZ Hz s 

and shows a tendency to decrease as load further increases. On the other hand, 
the 3D-rigid asperity curve shows a steady increase in ^ with P„ even for very 
heavy loads. The lack of agreement at these heavy loads is definitely due to the 
local elastic effect. 

Referring to Fig. (4.8), with the ratio h^/R changing from 10“^ to 10“^, the 
discrepancy between the 3D-rigid asperity results and the 2D-elastic asperity is 
greatly enlarged. In this regime of extremely thin film, the local elastic effect 
is so overwhelming that one cannot expect any validity of the 3D-rigid asperity 
perturbation analysis. 
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(3) Effect of G 


The effect of pressure-viscosity dependence is examined by varying G from 100 
to 500 for " 10*^, -0.5, b « 1/32, and 0.3. The reason 

for using a mild pressure-viscosity dependence is because recent studies in EHD 
traction has indicated that the effective pressure-viscosity coefficient within the 
Hertzian conjunction is only a small fraction of those measured under static equi- 
librium. For this reason, it is believed that values of G around 500 should not be 
unreasonable to represent the effective pressure-viscosity dependence in the con- 
junction. It can be readily seen in Fig. 4.9 that the magnitude of increases 

s 

with G for both 3D-rigid asperity and 2D-elastic asperity cases. However, the 

rigid-asperity model is slightly more influenced by G than the elastic asperity 

model. It is also interesting to note that the pressure perturbation A is 

s 

depending exponentially on the pressure-viscosity coefficient. This is somewhat 
expected since pressure is directly portionally to the viscosity. 

(4) Effect of 

With .003, hjj/R = 10‘^, G = 100, b » 1/32 ,Cj/fao = 0-3. Fig. (4.10) shows 

the relation between A and X^, the location of the center of the asperity. 

s o 

Qualitatively, the effects on A due to rigid and elastic asperity have the same 

s 

trend. The magnitude of Ag increases as X^ moving toward the center of the con- 
tact. This trend is certainly expected, since a higher pressure level inevitably 
would lead to a stronger asperity interaction, hence higher pressure perturbation, 

(5) Effect of b 

Fig, (4,11) shows the effect of b, asperity size, on A for P = 0.003, 

S nZ 

hj/R = lO"'^, G = 100, Xg® -0.5, c^/h^ 0.3. From the two curves shown, a similar 
trend is seen to exist In the relation between and the asperity size for both 
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the rigid model and the elastic model* In both cases, ^ increases with the asperity 

s 

size; however, the increase in ^ is much larger in the case of rigid asperity , 

s 

model than that in the elastic- asperity model. 

(6) Effect of 

With P„ - .003, h /R = lO"^, G = 100, X_= -0.5, b = 1/32, Fig. (4.12) shows 

Hz O 

the effect of asperity height, c^/h^, on Ag f or both 3D-rigid-asperity and 2D-elastic 
asperity models. It is readily seen that the two curves almost coincide with each 
other and the magnitude of for both cases increases with c^/h^. As c^/h^ 
becomes larger, the two curves begin to divert slightly, ^ 

(7) Effect of V 

As defined earlier, y Is the ratio of half of the asperity length to half of 

the asperity width, known as the ellipticity ratio of the asperity. As the 2D 

elastic asperity model only describes straight asperity ridges. Fig. (4.13) only 

shows the results for the 3D-rigid-asperity model. It is seen that the magnitude 

of ^ increases with y as expected. When y is less than 1, ^ increases sharply 
s s 

with a small increase in y . As y becomes larger the change of becomes much more 
gentle. In fact, for y > 5, the change of y is practically negligible, 

(8) Effect of Slide to Roll Ratio S 

Fig. (4.14) shows the effect of S on A for P = .003, G = 100, h /R = 10 

S rlZ o 

-0.5, b = 1/32 and c^/h^= 0,3. For y = 1,2 and 10, Ag decreases when S increases 

from -2.0 to 0.3, and then increases when S increases from 0.4 to 2,0. These trends 

show that between S = 0,3 and 0,4, the perturbed pressure A reaches a minimum. 

s 

This phenomenon can be explained by Fig. (4,15) in which the perturbed pressures 
$ around the asperity center are plotted for S = 0.3 and S = 0,4. In the case of 
S = 0.3, the value of $ is mostly negative for x/b' ^ 0, and positive for x/b* S: 0, 
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In the case of S = 0.4, the pressure exhibits an opposite trend. Thus, one would 

expect that between S = 0.3, and 0.4, the trend for $ begins to reverse itself. 

With P„ , h /R, Xo, b, c. /h and Y having the same values as those used for 

Fig* (4*14), Fig. (4.16) shows the effect of G on S. For G = 100 and 500, the 

qualitative trends of versus S are the same. However, for G = 100, the minimum 

^ is located between S = 0.3 and 0,4, whereas, for G = 500, the minimum ^ is 
s s 

shifted to a position between S = 0.1 and 0.2, Thus when the magnitude of G is 
increased, the miniimim is shifted toward the asperity center. This agrees very 
well with Lee and Cheng's [5] results in which the G value is equal to 3180, the 
perturbed pressure is negligibly small when S=0. 

The pressure profiles obtained from the smooth film theory and the perturbation 
theory are compared in Fig, (4.17). Again, it shows that in the case of pure rolling 
(S=0), the perturbed pressure can be neglected. For |s| =2, the perturbed pressure 
is relatively more important. 



4,5 CONCLUSIONS 


The perturbed pressure distribution around an ellipsoidal asperity tip within 
an EHD line contact can be calculated by solving the perturbed Reynolds equation 
based on the assumption of a rigid asperity. The magnitude of the perturbed pres- 
sure Ag, was found to be a function of the following dimensionless variables: The 
Hertzian pressure, ^Hz* smooth film center film thickness h^/R, the pressure 

viscosity parameter G, the position of the asperity center X^, the size of the 
asperity b, the height of the asperity c^/hQ, the ellipticity ratio y and the slide 
to roll ratio S. A was shown to increase with an increase of P , G, b, c-/h , Y> 

S rlZ X O 

or Xo* However, it decreases as h /R increases. The manner in which A varies with 
J o s 

S is dependent upon the pressure viscoisty parameter G. For a large G, A is at 

s 

its minimum when S approaches zero (pure rolling condition); it increases as the 

magnitude of S increases. For a small G, the value of S at which A reaches a 

s 

minimum shifts from zero to some small positive values. 

A comparison was made between the results obtained by the perturbation analysis 
based on an ellipsoidal asperity for Y = 10 and S=2, and those obtained by the 2D 
asperity analysis. The comparison indicated that the perturbation analysis which 
ignored the local elastic effect yielded a very good approximation on the magnitude 
of the pressure fluctuation for certain cases, provided that h^/R ^ 10 
P ^ 0.003, G 100, c-/h 0.3 and b =- 1/32, Beyond these limitations, the 
perturbation analysis would over-estimate the pressure fluctuation. 


/ 

/ 
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Figure 4-3 One-Dimensional Grids for Smooth Contact 
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Figure 4-4 Two-Dimensional Grids for Rough Contact 
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Figure 4-5 Grid Spacing 
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Fig. 4-6 The Effect of the Nominal EHD Center Film Thickness, 

h. /R, on the Double Amplitue of the Perturbed Pressure, 
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Fig. 4-9 The Effect of the Pressure Viscosity Parameter, 
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Fig. 4-11 The Effect of the Asperity Size, b = b'/b,on 
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Fig. 4-12 The Effect of the Asperity Height, c,/h , on ii.* J • 
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Fig. 4-14 The Effect of the Slide to Roll Ratio, S = 2(u^-U2)/(Uj_+U2>, 

on/Jt»5iM<‘’^«>*»ifor a Three-dimensional Rigid Asperity with 
Pressure Viscosity Parameter ,o< E ' = 100 
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Fig. 4-15 The Perturbed Pressure, $ , Around the Tip of a Three- 
dimensional Rigid Asperity 
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Fig. 4-16 The Effect of Slide to Roll Ratio, S = 2 (u- -U 2 )/ (u^+u^) 
5 •”*»» ^ Three-dimensional Rigid Asperity with 

Pressure Viscosity Parameter, c^E' = lOO and 500 





CHAPTER V 


THE EFFECT OF SINUSOIDAL WAVINESS ON THE PRESSURE 
FLUCTUATION WITHIN THE HERTZIAN CONTACT 

5.1 Introduction 

In the previous chapter, the emphasis has been placed on the pressure 
perturbation around a single asperity. Such analysis .with a single 
asperity is useful only in determining qualitatively the effects of 
asperity geometry, lubricant property, speed, and load upon the per- 
turbed pressure amplitude; but it is not suitable in making quantita- 
tive predictions of pressure fluctuations within these contacts. In 
order to simulate the effect of continuous transverse ridges on the 
pressure fluctuation in an elastohydrodynamic contact, one may assume 
these ridges can be represented by a series of sinusoidal waviness. 

In the case the surface roughness is located on the stationary surface 
only, the pressure and deformation profiles become time- independent, 
and one can modify the method described in \_ 22 ~\ to solve for the com- 
patible pressure and film thickness profiles. The analysis in this 
chapter is confined to seek the EHD performance at the inlet of a 
Hertzian contact with a sinusoidal roughness on the stationary surface. 

5.2 Mathematical Formulation 

The roughness model is assumed to be continuous transverse ridges 
which can be represented by a sine wave on the stationary side of the 
contact. Thus the asperity height 6, is a function of x only. There- 
fore one would obtain the relations 

I « c^/h^ sin(2n nX) (5.1) 
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where 


(5.2) 




0 


8 ■ 6/h0 


n 



or 


1 

4^ 


X - x/b 

h “ center film thickness of the smooth-film EHD contact 
Cj^ » maximum asperity height 


b* <B half of the asperity width 

b « half of the Hertzian contact width 

X » coordinate axis along the sliding direction 


The Reynolds equation and the elasticity equation are respectively 


dX ^ ^ V 3 j 




14. - + 6 


(5.3) 


(5.4) 


^Hz 




(5.5) 


where 


U 


p - p/Ph.. - »,/P, * - 

2E'R * ^ “ **X^**®* hj “ + 8» 


-1 




5 - 5/l>» 


G - aE' 
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Using the same numerical scheme stated in [22], P(X), and U are 

solved for a given set of non-dimensional variables, namely nominal 
center film thickness Hertzian pressure P^^^ssure viscosity 

parameter G, maximum asperity height Cj^/h^ and asperity width b*/b. 

5,3 Discussions of Results 

The effect of the pressure viscosity parameter G, the maximum asperity 

height c^/h^, the asperity width b*/b, the nominal center film thickness 

h ^R, and the Hertzian pressure magnitude of the perturbed 

pressure, , (which is defined as the difference between the maximum 
s 

and the minimum pressure ripples deviated from the nominal smooth-film 
pressure profile), and the actual pressure distribution will be discussed 
in the following sections. 

5.3.1 The Effect of Pressure Viscosity Parameter G 

For Pjj^= 0.003, h^R = lO"^, = 0.3 and n = 2(b'/b = 1/8). Fig. 

(5.1) shows the effect of G on A . Similar to the results obtained for 

s 

the single asperity ridge in an EHD contact, as discussed in Chapter IV, 

the magnitude of ^ increases with G for this case with a waviness sur- 

s 

face profile. With the same set of values for P^^’ ^l^^o 

b*/b, the effect of G on the pressure profile is shown in Fig. (5.2). 

As expected, the pressure fluctuation from the nominally smooth-film 

pressure profile is more pronounced ^en the value of G is larger and 

within the Hertzian contact zone. Since G represents a significant 

parameter, the effects of other parameters on A are compared in the 

s 

subsequent sections for G = 100 and G = 1,000. An attempt was made for 
cases with G greater than 1,000, some of the solutions failed to converge. 


81 



5,3,2 The Effect of c-/h 
1 o 


The asperity interaction will be more severe when the amplitude of 
the asperity Is larger. This phenomenon has been shown earlier In the 
single asperity-ridge analysis and is exhibited again here in Fig, (5.3) 
for the case of wavy surface roughness with P„ = 0.003, h /R = 10 and 
n ** 2(b'/b * 1/8), For both cases with G = 100 and 1,000, the magnitude 
of the perturbed pressure increases with the ratio Cj^/h^. This phenomenon 
can also be observed in Fig. (5.4) with G = 100 and Fig, (5.5) with 
G ® 1,000 in which the perturbed pressure profiles are shown. 

5,3.3 The Effect of n 

For P„ = .003, h /R = lO"^ and c,/h =0,3, Fig. (5.6) shows the effect 
tiz O 1 o 

of n on the magnitude of the perturbed pressure ^ , for both G = 100 and 

s 

1,000. It seems that the two curves shown here do not have the same 
qualitative trend. However, it is believed that the further increase of 

n for the case with G = 1,000 will decrease the magnitude of ^ due to 

s 

the elastic effect, so that the shape of the curve in this case will be 
similar to that with G = 100. 

Let's recall that n = 1/4 xb/b'. Ihus increasing the value of n 
implies the closer the distance between the center of each individual 
asperity ridge and the center of the EHO contact. Therefore the effect 
of n actually consists of the effects of elastic deformation, the ratio 
b'/b and the distance between the individual asperity center and the 
contact center. 

The perturbed pressure profiles for G = 100 and 1,000, due to the 
effect of n, are shown in Figs. (5,7) and (5.8) respectively. The 
results obtained for these two cases are the same qualitatively. For 
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G *= 1,000 in which the elastic effect is more severe, the magnitude of 
the pressure ripple deviated from the nominally smooth-film pressure 
profile is much larger than those with G = 100. 

5.3.4 /R Effect 

The effect of h /R on the magnitude of the perturbed pressure ^ , for 
o . . s 

both 6 100 and 1,000 are shown In Fig. (5.9). In these exampies, the 

values of P„ , c,/h , n and b'/b are 0.003, 0,3, 2 and 1/8 respectively, 

HZ i O 

It is observed that the magnitude of ^ increases with the ratio h /R. 

,s o ^ 

This increase is larger when the magnitude of G is larger. 

In Chapter IV, it has been shown that the magnitude of ^ decreases 

s 

with the increase of h^/R when there is a single 3D rigid asperity within 
the contact zone. This opposing trend between the case of a single 3D 
rigid asperity and continuous elastic asperities can be explained in the 
following manner: 

In ^he case of sinusoidal elastic asperities, the pressure amplitude 

is reduced by the local elastic deformation of the asperity. If c-/h 

I o 

is held constant, the reduction in pressure due to local elastic defor- 
mation becomes much greater as h /R decreases, because at a smaller h /R, 

o o 

it requires only a very small pressure amplitude to flatten out the 
asperity. This phenomenon of elastic effect is best illustrated by 
Fig. (5.10) which shows the perturbed pressure and surface profiles for 
h^/Rr = 10 5 X 10 ^ and 10 ^ respectively, for the sinusoidal elastic 

asperities with G = 100. The shapes of these perturbed pressure and 
surface profiles are consistent with one another. However,, the smaller 
the h^/R, the more the asperity being flattened out. For the case of a 
3D rigid asperity, the effect of local elastic deformation is ignored, 
and thus the opposing trend is found. 
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nil III 


II I 


The perturbed pressure and surface profiles for the sinusoidal eleistic 

asperities with G =* 1,000 are shown in Fig, (5,11). The results for 

both G <= 100 and 1,000 are shown to have the same characteristics and 

the magnitude of A is larger with G = 1,000 as expected, 
s 

5,3,5 P„ Effect 

ijz 

For h /R = lO”^, c,/h = 0.3, n = 2 and b'/b •* 1/8, the results of 
o i o 

the magnitude of the perturbed pressure A , versus the nominal value of 

s 

the maximum Hertzian pressure P , are shown in Fig, (5,12). It is 

£lZ 

readily seen that, for both cases with G = 100 and 1,000, the increase 
in P„ results in a decrease in A . This phenomenon is due to the 
reduction of the pressure amplitude caused by the local elastic defor- 
mation of the asperity. The larger the P , the easier the asperity 
being flattened, and thus the smaller the 

The perturbed pressure profiles are shown in Figs* (5.13) and (5.14) 
respectively. The results for G = 100 and 1,000 are found to have the 

same trend qualitatively. As expected, the magnitude of A is larger 

s 

in the case of G = 1,000. In addition the perturbed surface profiles 

for the case of G = 100 as shown in Fig, (5,13) illustrates the effect 

of the local elastic deformation of the asperity as explained previously. 

It tells that the larger the P , the more severe the asperity being 

Hz 

deformed. 
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5.4 CONCLUSIONS 


In the case of the sinqple sliding of a smooth surface against a 
stationary rough one, the magnitude of the pressure deviation from the 
nominally smooth-film profile and the perturbed pressure profile in the 
inlet of an elastohydrodynamic contact, can be determined quantitatively 
When the undeformed rough surface profile is given. 

In the examples given in this chapter, the undeformed rough surface 
profile is simulated by continuous transverse ridges represented by a 
series of sinusoidal waviness. The effect of the following non- 
dimensional parameters on the magnitude of the pressure fluctuation L 

s 

and the perturbed pressure profile are obtained: the pressure viscosity 
parameter G, the maximum height of the asperity Cj^/h^, the number of 
wave cycles within the contact zone, n, (which In turn determines the 
width of the asperity, b'/b), the nominal smooth-film center film thick- 
ness h^R and the Herti^lan pressure 

The effects of G and magnitude of for this case with 

a wavlness surface profile have the same characteristics compared to 
those obtained for the single asperity ridge In an EHD contact as pre- 
sented In Chapter jCV. Namely, the magnitude of ^ Increases when the 

s 

magnitude of G or Cj^/hp Increases, 

The results of versus P„ and h /R are found to be consistent for 

■ 0 H.Z 

both G >= 100 and 1,000, The magnitude of ^ decreases as the magnitude 

s 

of P„ or h /R becomes larger. These results for the waviness surface 
Hz o 

profile have the opposing trend when compared to those obtained by the 
3D rigid asperity analysis. The reason for the opposing trend Is due to 
the pressure reduction caused by the local elastic deformation of the 
asperity for the wavlness surface profile, whereas the effect of local 
elastic deformation Is Ignored for a 3D rigid asperity. 
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The effect of the number of wave cycles within the contact zone on 

the magnitude of the perturbed pressure t Is found to be very compll- 

s 

cated. The local elastic effect, the ratio b'/b, and the distance 

between each Individual asperity center- and the contact center are 

Important factors affecting the magnitude of A caused by changing the 

s 

magnitude of n. 
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the Perturbed Pressure Profiles for Pressure Viscosity 
Parameter, E ' = 100 
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THE SMOOTH- FILM PRESSURE 
AMD SURFACE PROFILES 
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Fig, 5-14 The Effect of the Normalized Hertzian Pressure, P = p , 
on the Perturbed Pressure Profiles for ot E' = 1000“^ ^ 
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APPENDIX A 


Justification of M as a Random Quantity of Zero (or Negligible) Variance 

When ^ is assumed to be zero, Eq. (2.6) becomes 
.0 1 


(A.l) 


where M is defined by Eqs. (2,7) and (2.10). Hence M is a constant of x. However, 
this constant can be a stochastic quantity in the time domain. 


EHD Contact 

Re-arrange Eq. (2.10) to obtain 

u, + u 


1 2La = iL IL 

12^, Bx 3 ■ 


2 h__ " ^2 

2 3 2 > ^ 




(A. 2) 


At X * - 00 , the asperity effect is negligible so that q = p = 0. However, the 
pressure p, must reach a significant fraction of the Hertzian maximum at x = - 1, 
such that e~^^ «1, i.e, Therefore q at x = - 1, can be assumed to be a 

stochastic variable with extremely small variance. Integrate Eq, (A. 2) from 


xe-ooto- 1, and obtain 
-1 

5TS ■ J 3 




“l +"2 


-1 




dx + 


Ui - U2 


-1 


J ^ 


6, - 6, 


dx 


(A. 3) 


As M is a constant of x, Eq. (A.3) can be re>wrltten as 

-1 


1 1 ^ ^ f ^ ^ ^2 . ^ “2 r 


u„ h + t-i 


dx 


M 


-1 


(A.4) 




dx 


When there are enough asperities within the contact zone, the integrals in Eq. (A.4) 
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are independent of the precise arrangement of the roughness and also independent 
on time. Therefore M can be assumed to be a stochastic quantity with zero (or 
negligible) variance. 


Rigid Rollers 

Eq. (2,7) can be re-written as 

^ -"^2 h , - ^2 

12y, dx , 3 " 2 3 2 



Integrate Eq. (A. 5) with the boundary conditions: 


p = 0 


at X = - cx> 



* 

at X = X 


(A. 5) 


(A. 6) 


one obtains 


* 

rtX h 


6i + 62 


2 P 

dx + 


* 

X h + 6 - 


dx 


M = 


r 


1 



dx 


(A. 7) 


Even though x is a random quantity, its deviation from the mean is expected to 
be of the same order of the wavelength of the asperity, and is small compared 
with the size of the bearings. The integrals in Eq, (A, 7) are only slightly 
affected by x , and therefore can be considered as stochastic quantities with 
negligible variances. Thus, M is also a stochastic quantity with negligible variance. 
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Inf Inttely-Wide Sltder 


Integrate Eq. (A. 5) with the boundary conditions: 


p *= 0 at X «= 0 and L 


(A.8) 


one obtains 


M = 



h - 6i + #2 



o 



dx 


(A. 9) 


Following the above argtraient, M is also a stochastic quantity with zero (or 
negligible) variance. 
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APPENDIX B 


FORTRAN IV LISTINGS OF COMPUTER PROGRAMS 

1. PROGRAM ROLSLIP: This program which is used in Chapters II and III, 

is to calculate the integrated pressure at the inlet of an EHD 
contact with random surface roughness by the stochastic approach. 

2. PROGRAM RIGR0L7; This program which is used in Chapter II, is to 
compute the load of a rigid roller bearing with random surface 
roughness by stochastic approach. 

3. PROGRAM SLIDERS: This program which is used in Chapter II, is to 

compute the load of an infinitely-wide slider bearing with random 
surface roughness by stochastic approach. 

4. PROGRAM DAP0L2: This program which is used in Chapters II and III, 

is to generate data for functions G 2 , G^ and G^ for the correspond- 
ing These data are input to programs ROLSLIP, RIGR0L7 and 

SLIDERS. The data for functions G 2 , G^ and G^ in the case when 
the asperity height distribution functions are in the form of 
polynominal function as well as sinusoidal function are tabulated 
respectively. 

5. PROGRAM ROLWAVl: This program which is used in Chapter III, is to 
compute the integrated pressure at the inlet of an EHD contact with 
waviness surface roughness by the deterministic approach. 

6. PROGBIAM ASPERITY: This program which is used in Chapter IV, is 

to compute the perturbed pressure distribution due to a single 
3D- Rigid Asperity. 
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PROGRAM R0L5LiH( INPUT, OOTFOT) 

C THIS program is to CALCULATE Th£ INTEGRATED PRESSURE AT THE INLET OF AN EKQ 
C CONTACT C5Y THF STOCHASTIC APPROACH 
C THIS PROGRAM CONSISTS OF FIVE SUBROUTINES 

C 1 - SUBPOUIINE GRUOIN IS TC CALCULATE THE INLET INTEGRATED PRESSURE OF AN 
C EhD CONTACT BY SMOOTH FILM CRUBIN THEORY 

C 2 - SUBROUTINE CONS IS TO CALCULATE THE CONSTANT C FOR VARIOUS DSIG DATA 

C 3 - subroutine FUNC is TO EVaLUaTe FUNCTIONS FI AND F2 FOR DIFFERENT S/h 

C 4 - SUBROUTINE TABLE IS TO SET UP FUNCTIONS FI AND F2 IN THE FORM OF ARRAYS 

C 5 - SUBROUTINE INT IS AN INTEGRATION ROUTINE TO CALCULATE THE INTEGRATED 

C PRESSURES WLOADl* WLOAD2, AND WL0AD3 
C 

C THE’ PRESSURE IS INTEGRARED FRCM X=-5 TO X=-l. BECAUSE OF THE NONLINEAR 
C PROPERTY OF THE INTEGRAND NEAR A=-l, THE INTERVAL OF INTEGRATION IS DIVIDED 

C INTO T»<C PARTS, I.E. FROM x=-5 TO X=-2, AND FROM X=-2 TO X*-i. 

C 

C DATA CARD UO. 1 
C KF = total no. OF GRIDS 
C r.M = GRIO NO. AT X=-2 
C NSIG = NO. OF OSiG(I) DATA 
C NHO = number OF HO DATA 
C DATA CARD NO. 2 
C w = LOAD PER unit WIDTH 
C SLIP = SLIPPAGE COEFFICIENT 

C PI = 3.1‘t1593 

C DATA CARD NO. 3 

0 OSIG(l) = DATA FOR SIGMA/HO 

C DATA CARD NO. 4 
C UhO(I) = VARIOUS HO DATA 
C DATA CARD NO. 5 

C THESE data ARE THE OUTPUT FROM PROGRAM SUCHAS PROGRAM DATP0L2 WHICH EVALUATE 

C The integral of t distribution FUNCTI0N*HSS>/ ll. + SH*HSS)**3 AND 
C (DISTRJBUTICN function/ ( l.+SH<»hSS)**3 

C This DISTRISUTIOi'i FUNCTION CAN BE POLYNOMIAL , GAUSSIAN, SINUSOIDAL OR ANY 
C other FUNCTION 

C OSH(I) = THE ABSCISSA RANGING FROM 0 TO 0.333 
C OV(I) = THE NONDlMF.NSIONAL DSH(l) 

C DGh(I) = THE INTcGRAL OF 35 /9b" ( 1 . -HSS**2/9 ) **3/ ( 1 ♦SH*HSS ) **3*HSS 
C IN THE CASE OF POLTNCMIAL DISTRIBUTION 

C OG2(I) = THE INTEGRAL OF 5 S /96" ( 1 . -HSS**2/9 ) *"3/ ( 1 ♦SH*HSS ) *»3 
C IN THE CASE OF POLYNOMIAL DISTRIBUTION 

C 065(1) = DG2(I)/uG‘»(l) 

C 

C THE OUTPUT OF THIS PROGRAM WGRUB, WLOADl, WL0AD2, WL0AD3, Wl, W2, W3 ARE ALL 
C INTEGRATED PRESSURE AT THE INLET 
C 
C 

COMMON WLCADi (I**) ,wLuAD2 (14) ,y.L0AD3 (14) , DS 1 6 ( 14 ) , T 1 ( , 301 ) » 

1 T2(1a»301 ) ,XX (301 ) .KK.hF »H0,W,P1,S,X,H 1,F2,C, SLIP, 

2 G(301 ) ,F,WGRLil1(1a) ,OG‘t(2Cl) ,0G2(20l) 

DIMENSION Q1 ( IH, 301 ) ,02 ( 14,301 ) .03 (14,301 ) ,DSH(201 ) ,DV (201) , 

C DG5(20l), DH0U4), WKii*), W2(14), W3(14) 

C 
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^ cvj ro ^ tn ^ 


PRINT i 

READ £»KF»KM*NSIG»NH0 

READ 3.h»SLIP*PI 

READ Jt (DSIG(I) »I=l»NSIG) 

READ (DHO(l) ♦I=l«NhO) 

DO 100 I=l»201 

100 read 5»DSh (I ) ’OV ( 1 ) «GG4 ( I ) »OGSn ) «062{I) 
rPRlNT 51 

'PRINT 2, KF, KM, NHO 

print 52 

PRINT *»*W»SLIP»PI 
PRINT S3 

PRINT A. (DSIG(I) ♦I=1»NSIG) 

PRINT 5a 

PRINT (OHOU) *I=1*NH0) 

XX(i) = -5, 

DO 98 l=2»KM 

98 XX(I) = XX(I-l) ♦ O.Olb 

KMM = KM ♦ 1 

DO 99 I=KMM,KF 

99 XXU) = XX(I-l) ♦ 0.01 

PRINT 56* 

DO 97 I=i»NHO 
IGRuB = I 
HO = DHO<I) 

CALL GRuBIN(IGRuB) 

97 PRINT 7,I,DHO(l>,W6RUea) 

DO 101 LL =1*NSIG 
PRINT b5.DSIG(LL) 

PRINT 6 
CALL CONS(LL) 

DO 103 I=l»NHO 
KK = I 

hO = DM0 (KK) 

S = DSIG(LL) « HO 
CALL TABLE 
CALL I NT 

WKI) = XLOADKI),/ WGRUE(I) 

W2(I) = kiL0AD2(I) / WGROc(I) 

W3(I) = WL0A03(I) / WGKOBd) 

PRINT 7,I»DH0(1) »WL0A01 (I) »WL0A02(I) »WL0AD3(I) »tal (I) *W2(I) ♦W3(I) « 
7 wGPUe(I) 

103 CONTINUE 
ici continue 
C 

FORMATdHl ) 

FORMAT (2015) 

FORMAT (6F10, 8) 

FORMAT (9(2X»E13.6) ) 

FORMAT (F10.6»Fi0.2»3F20.10) 

FORMAT (/3X»»Ia»6X«*0H0»»9X»»wLUAo1**9X«*WL0AO2*»9X**WL0Ao3**13X» 

6 *wl«*i3X»*W2*»13X»«*l3*,lCX,*WoRUB*) 
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7 F0kMATtI5»8(2X»E13.6)) 

51 F0KMAT(/3X»*KF*. JX»*KM*,X»*NS1(j*»2a»*NH0*) 

52 FORMAT (/7Xt*t»l*»llX»*SLlP«»13X»*PI*) 

53 format 1/6X»*DS16(I)«) 

54 FORMAT (/6X»*dhO(I)*) 

55 FORMAT (///6X**uSIG=**F10,e/) 

56 FORMAT (/3X»*I»»6X»*DH0*»9X ,*WGRUB*) 

• STOP 

ELfvD 

SUBRObllNE GRod 1 N ( 1 OKUtt ) 

COMMON WLCADl (14) *'.vL0AD2 ( 14) ♦wLOA03(l‘«) »OSIG (l^) »Tl <14,301 ) * 

1 T2(14»301),XX(3C‘l)»Kr;»KF»KM,NF,HO,W,Pl,S,X,Fl»F2,C*SLIP* 

2 G(301) ,F,ii/GRUb(14) »DG4(2C1) ,062(201) 
dimension F3(j01) 

c 

DO 400 J=1,KF 
X = XX(J) 

H = H0 + 4»*lrvPI*(ABb(X)*syRT ( AdS ( < “X ) **2-l .0) )-ALOG(AbS(X) +SC)RT ( A6s 
1 ( (-X)**2-1.0) ) ) ) 

AOO F3(J) = (H-HO) / H**3 
WGRUB(iGRUB) = 0, 

KMK = KM - 2 
DC «01 J =1»KMK,2 

401 wGRUe(lGRUB)=WGRUB(IGRU8) ♦ 0 . 0 1 5/3 . * ( F3 ( J ) ♦4,*F3 ( J + 1 ) ♦F3 ( J>»2) > 
KFF = NF - 2 

DO *+02 J = KM,KFF,2 

402 wGRUB(iGRUB)=WGRUB(IGRUb)+ O.Ol /3 .* (F3 < J) +4 ,«F3 ( J+ 1 ) +F3 ( J^2 ) ) 
RETURN 

END 

bUBROU t INl CJiOlLLJ ■ “ 

COMMON wLCADl (14) «WL0A02(14) ,wL0AD3(14) ,DSIG(14) ,Tl (14,301)* 

1 T2(1-4»301) ,X^(301) ,KK,KF*KM,NF ,H0,W,PI,S,X,F1,F2,C,SLIP* 

2 6(301) »F,WGRUd(l4) ,DG4(ZC1) ,DG2(201) 

C 

V = 1* ♦ USIG(LL) / 0.001665 
MV = V 
VM = MV 

G2 = (l.'G2(MV*-1) - DG2(MV)) * (V-VM) ♦ DG2(MV) 

64 = (uG4(MV+l) - DG4(hV)) * (V-VM) ♦ UG4(MV) 

C = G4 / 62 
RETURN 

E>iO 

SUbROUilNE TABLE 

COMMON WLOADl (l4) »wL0A02(l‘+) ,WLOAD3(14) , DS IG ( 14 ) , T 1 ( 14 , 30 1 ) * 

1 T2(14,301) ,XX (301) ,KK,KF,KM,NF.Ha,W,PI,S,X,Fl,F2,C,SLIP* 

2 G(301) »F»WGRUo(l4 ) ,UGh(2C1) ,DG2(201) 

DO 201 J=1»KF 

X = XX ( J) 

call func 

T1(KK,J) = FI 
T2(KK,J) = F2 
201 CONTINUE 

return 

END 
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c 


c 


c 


c 


301 


302 


C 


"SUbf^OUTlNfc. 

COMMON iVLCADl «f.^0AD2(i4) f hL0A03 (l4),DSIG(14)*Tl(l<t*301)f 

1 T2(1A»301 > .XX (301 ) .KK.KF .KM.NI- .HO.W. PI.S.X.Fl .F2»C»SLIP» 

2 0(301) .F.WGRUb(IA) .DGA( 2C1) .DG2(20l) 


H=H0*A.«W/PI« (ABS(X)*SQRT (ABS( (-X) **2-1.0 ) )-ALOG(A0S(X) *SQRT(AiS 
1 ( (-X)<^*2-1.0) ) ) ) 

SH = S / H 

V = 1. ♦ S/H/0. 001665 
MV = V 
VM = MV 

G2 = (l»G2(MV + 1) - DG2(MV)) * (V-VM) ♦ 0G2(MV) 

GA = (UGA(MV+l) - DGA(MV)) * (V-VM) ♦ DGA(MV) 

FI = (ri-HO) /H*<>3*G2 

F2 = (GA - C*-G<i) / h<f*3*S 

RETURN 

END 

SUdROUiINt JnT 

COMMON WLCADl (I*.) . WL0AD2 ( lA) .wL0A03(lA) »DSIG ( 1 A ) » T 1 ( lA . 301 ) . 

1 T2(lA.301).XX(301).KK.KF.NM.NF.hO.W.PI.S.X.Fl.F2.C»SLIP. 

2 G(301) .F.WGRUnd^*) .UGa(2C1) .DG2(201) 
dimension PI ( 1*.) ,P2 (1a) 


PKKK) = 0. 
P2(KK) = 0. 

KMK = KM-2 
00 301 J=1.KMK.2 


PKKK) = FI(KK) + 0. 015/3. *(T1 (KK.J)+A.*T1(KK.J + 1)-.T1(KK»J*2) ) 

P2(KK) = P2(KK) ♦ 0. 015/3. *(T2(KK.J)*A.*T2(KK.J*1)^T2(KK.J-.2)) 

CONTINUE 

KFF = kF-2 

DO 302 J=KM,KFF,2 


PKKK) = PKKK) + 0.01/3.« (TK^K,J) +A.*T1(KK, J*l) +TKKK,J*2) ) 
P2(KK) = P2(KK) ♦ 0 .01/3.* (T2 (KK, J) ♦A.*T2 (KK. J+1 ) +T2 (KK» J*2) ) 
CONTINUE 

WLOaDKKK) = Pi(KK) - P2 (kK)*SLIP 
WLOAD21KK) = Pl(KK) ♦ P2(KK)*SLlP 
KLOAD3 (KK ) = PKkK) 


RETURN 

END 

OOOO&OCOOOOui^OOOOOOOOO End of PECORU 
301 201 9 6 

0.00003 0.1 3.14159265 

0.05 0.1 0.1b 0.2 0.25 0,3 0.31 0.32 

0.33 

O.OuOOl 0.00002 0,00003 ti.0o0u5 0.00007 0.00009 

C DATA FOR DSH(I). Dv(I). DC--t(l)* ag2 ( I ) . Dg5 ( I ) AWE OMITTED 

ooooooooooooooooooooooen'D of Information 


PROGRAM RIGROLMINPuTtOUTPOT) 


C 


C THIS program is to COMPUTE TP£ LOAU OF A RIGID ROLLER BEARING BT STOCHASTIC 
C APPROACH 
C 

C THIS program CONSISTS OF A SUBROUTINES 

C I - SUBROUTINE TABLE IS TO SET UP A PRESSURE ARRAY CORRESPONDING TO THE GRID 
C POSITIONS 

C 2 - SUBROUTI/JE CONS IS TO CALCULATE THE CONSTANT C 

C 3 - subroutine FUNC is to EVALUAfE FUNCTIONS Fit F2t AND F 

C *» - SUBROUTINE INT IS A SIMPSON imEGRAT ION ROUTINE TO COMPUTE THE LOAD W(KK) 
C THE OUTPUT OF THIS PROGRAM IS w(I) AND THE CORRESPONDING DSIGJI) 

C 

C DATA CARD HO. 1 

C M = NO. OF uniform GRIDS FROM X=-l TO X=0 

C NF s NO. OF UNIFORM GRIDS FROM A=0 TO X=XA 

C KO = M ♦ 1 
C KF = NF ♦ I 

C ND = NO. OF D5IG(I) DATA 
C data card i'lO. 2 

C XA(i>.Xa(2T initial values ASSIGNED TO XA(I)t ASSUMING THAT THE PRESSURE AT 
c This location is o. 

c HO = center film Thickness 

C SLIP = SLIPPAGE COEFFICIENT 

C PI = 3.1H1593 


C DATA CARD NO. 3 
C OSIG(l) = SIGMA / HO t DATA 
C DATA CARD ImO. A 

c These daia are the output from program suchas program oatpolz which evaluate 
C THE integral OF ( DISTfilbUTIUN FuNCT ION*hSS) / ( 1 . tSH^nSS) **3 AND 
C (DISTRIBUTION P UNCT I ON/ ( 1 . ♦SH*HSS ) **3 

C This DISTRIBUTION FUNCTION CAN BE POLYNOMIAL tGAUSSIANt SINUSOIDAL OR ANY 
C OTHER KUNCTION 

C OSH(l) = THE ABSLISSA RANGING FROM 0 TO 0.333 
C UV(I> = The NONDIMENSIONAL OSHU) 

C UG-Cl) = THE integral OF 3£ /9fa* ( 1 . -hSS**2/9 ) **3/ ( 1 tSH^bSS ) **3*HSS 
C IN THE CASE OF POLYNOMIAL DISTRIBUTION 

C DG2(I) = THE integral OF 3 5/ 96^^ 1 1 . -HSS**2/9 ) **3/ ( 1 *SH*HSS ) **3 
C IN THE CASE OF POLYNOMIAL DISTRIBUTION 

C DG5(I) = DG2(l)/OGA(I) 

C 

COMMON XX (121 ) »P (20» 121 ) tCSlGC 10) tW (10) .HOtHAtKOtKFtK.KKtSOtSFtXt 
C F.Fl,P2tDELA.Pl.C.CGA (201 ) tUG2 (201 ) .SLIP 
DIMENSION XA (2i ) tDSH (201 ) tUV (201 ) .UGS(20l) tXB U 0 ) t IK ( 10) 


C 


95 


print 1 

READ 2.NI tNFtKO.KFtND 

READ A.XA (1 ) tXA (2) thOtSLIPtPl 

READ A, (OSIG(l) tl=l.NO) 

DO 95 1=1.201 

read P.DSH( 1) tOV ( I ) .06*4 (I ) t065( I) .DG2(I) 
print 51 

PRINT 2.NI .NF.KOtKF.NO 
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u» t\) ►- r> 


PRINT t>2 

PRINT V*XA(l),AA<2>*H0tSLIP*Pl 

print o 2 

PRINT (DSIG(I) ♦I=i*ND) 

SO = l./NI 
DO 99 0=1 tKO 

99 XX(J) = -1, ♦ (J-1)*S0 
DO 100 II«1,ND 

KK = II 

print o,ii»DSiO(ii) 

XA<1) = 0.01 
XA(2) = 0.02 
DO 101 I=l»2 
K = I 

HA = 1. ♦ 0.5«XA<I)*»2/hC 
SF = XA(I) / nF 
PRINT bA 

PRINT J*I*XA(I) ttiAfSOfSP 

call cons 
CALL table 

101 CONTINUE 

XA(3) = (XA(1)*P(2*KF) - X A (2) *P ( 1 *KF ) ) / (P(2*KF) - PdtKFl) 
DO 102 I=3»20 
K = I 

HA = 1. ♦ 0.5*XA(I)*«2/HC 
SF = XA(I) / NF 
PRINT b4 

PRINT 3*ItXA(I)»HA.S0»SF 
CALL CONS 
CALL table 

DIFF = P(I*KF) / P(I»KO) 

IF (ABb(DIFF) .LE. O.Oul ) GO TO 103 
PROD = P(I»KF)"P(I-1»KF) 

IF (PROD .GE. 0.) GO TO 104 

XA(I*1)= (XA(I-1)*P(I»KF)-XA(I)*P(I-1*KF) )/(P(I*KF)-P(I-l»KF>) 
GO TO i02 

104 CONTINUE 
P(I-rl.KF) = P(I-2*KF) 

XA(I-l) = XA(I-2) 

XA(I+1)= (XA(I-1)*P(I»KF)-XA(I)*P(I-1*KF) )/(P(I*KF)-P(I-1*kF>) 

102 CONTINUE 

103 CONTINUE 
CALL INT 
X8(II) = XA(K) 

IK(II) = K 

100 CONTINUE 
PRINT 7 

DO 105 II=1»ND 

105 PRINT ti»II.DSIU(II),W(II)»X3(II)»lK(II) 

FCRMAT(lHl) 

FORMAT (2015) 

F0RMAT(I5*&(2X»E13.6>) 

FORMAT(eFlO.B) 
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5 format (F10.6^Fi0.ii«'jFa0.1C) 

6 FCt'VAT (//////X»«I I=«» IH, 1CX,<>L)SIG=*»F8.4//) 

7 FORMAT (////3X.*II*» 10X»*CS i6X»-»kv»»19X»^XA'>* 12A»*K«/) 

8 FCkmaT(I5»3(3X<E17.10) »3X»I2) 

9 FOf^MAT (9 (2x.El3.fe) ) 

51 FOKMAT 1/3X.<»NI». JX»«NF<^»3X . *K0* ♦ 3X . *KF* ♦ 3X » «ND» ) 

52 FCkMaT (/6X ♦'-i'XA < 1 ) *♦ 10X»“XA ( 2 > . 1 3X ♦ «HO* . 1 i X » *SL IP* ♦ 13X ♦ *P I* ) 

53 FOkMAT l/feX.oC'SiGd)*) 

54 format (/4X» i»r>»6X.*X A ( I ) « . 1 3X » *H A* . 1 3X ♦ * 50 * ♦ 1 3X » *SF* ) 

C 

RETURN 

END 

SuBwOL'IlNh TABlL 
C 

COMMON XX(121) *P(20.121) .CfelGllO) »W(10) »HO*HA»KO»kF.K»KK»SO*SF*X« 
C F »F1 2.DELA.P1 »C. DGa (201) .DoA (201 ) ♦ SLIP 
dimension E(12i) 

C 

PRINT 21 
1 = 1 

X = XX (J) 

CALL FuNC 
t ( J) = F 
P(K,J) = 0. 

DO 201 J=c»KO 
X = XX (J) 

CALL F.UNC 
E ( J) = F 

X = (XX ( j-i ) .XA ( J) )/2. 

CALL FUNC 

P(K,J) = P(K«J-1) ♦ SC/2 . / 3 . * ( £ ( J- 1 ) *4 ,*F *E ( J ) ) 

201 CONTINUE 
KOI = KO + 1 

DO 202 J=K01.kF 
XX(J) = XX(J-l) ♦ SF 
X = XX (J) 

CALL FuNC 
E(J) = F 

X = (XX(J-l)+xx(J) )/2. 

CALL FUNC 

P(K,J) = P(K.J-l) ♦ SF/2./3.*(£( J-1) +4.*F*E(J) ) 

202 CONTINUE 

DO 203 ‘J=1»K0.4 

203 PRINT 22.J.XX(J) »P(K.J) 

DO 204 J=KC1.KF»2 

204 PRjNT a2 ♦ j ♦ XX ( J) »P (k * J) 

C 

21 FORMAT (//3X »*j* » 1 OX ♦ *XX« . 1 8X » *P*/) 

22 FORMAT (I5.2(3X»E17. 10) ) 

RETURN 

END 


/ 
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SUtsROUfINt CONb 

COMMON XX (121) »P(Z0»l£l) tCSIGUO) «w(lO) ♦HO»HA»KO*KF»K*KK*SO»SFtXt 
C F»F1 ♦F2»DELA*i'IfC»U6‘* (<-01 ) »DGi: (201 ) ^SLIP 

V = 1. ♦ DSIG(NK) /HA/0. 001665 
MV = V 

VM s MV 

G2 = (0G2(MV>1)-DG2(HV) )«(V-VM) ♦ DG2<MV) 

Ga = (0G4(MV*1)-DGA(MV) )« (V-VM) ♦ DG4(MV) 

C = G4/G2 
PRINT 4l,K.K»K»C 

41 FORMAT (/3X»*KK<^»^X,*K*»1CX»*C»/2I5»3X«E17. 10/) 

RETURN 

EM) 

SL'HROUl INt FUNR ' ~ 

COMMON XX(121) ♦P(20»121) tCSIGdO) .i»(10) ♦H0»HA»K0.KF*K»KK»S0»SF,X» 
C F*Fl«f-2.DFLA»RltC*OG4(^:Oi) *062 (20 1 )» SLIP 
C 

H = 1. ♦ 0.5*X<»«2/HO 

V = 1. ♦ CSIG(NK)/H /O. 001665 
MV = V 

VM s MV 

G2 = (UG2(MV+1)-DG2(MV) )«(V-VM) ♦ DG2(MV) 

G^ = (UG4(MV*l)-OG‘t(Mv) )«(V-VM) ♦ 0G4(MV> 

Fl = (h-HA) /H**3«G2 
F2 = (oi,-C*G2)/h**3 
F = Fl ♦ CSIG(NK) * F2 « SLIP 

return 



SUtiROUlINt InT 

COMMON XX(12l)»F(20*12l)»CSIG(i0)*w(10)*HO»hA»KO»KF»K»KK»SO»SF#Xt 
C F»F1»)^2»DELA«P1»C»0G^{2C1) tUG2 (2ol ) »SLIP 
C 

W(KK) = 0. 

KC2 = KC>2 

DO 30 1 J - i , K02 , 2 

WtNK) = W(KK) ♦ SO/3.* (R (N» J) ♦ 4.»P(KtJ+l) + P(K.J+2)) 

301 CONTINUE 

KFF = NF - 2 
00 302 J =K0*KFF*2 

*lt(KK)'= W(KK) ♦ Sr/3,«(P(t<» J) ♦ M.*P(KtJ*l) ♦ P(K»J*2)) 

302 CONTINUE 
PRINT jl»NK.W(NK) 

C 

31 format (////X,*Kk=*, 13. i0X»*W=»»E17.10//) 

RETURN 

END — r 

OOOOuOOOOOCOuOOOOOOOUO END uF htCOHO 
Hi) 10 41 51 10 

0,01 0.02 0.0005 -C.l 3.14159265 

0,u6 0.1 0.14 C.ia 0.22 0,26 

0.32 0.33 

C data for DSH(I)< DV(I). UG‘«<I)» 4G2 ( I ) » DG5(I) ARE 
OOOOOOOOOOOOOOOOOOOOuO EM) Ul INFORi'iATION 


0.3 0.31 

OMITTED 
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program SLIOER3(INPuTfOUlFUT) 

C ■ j 

C THIS program is to compute ThE LOAL) OF AN INF INI TEL Y-WIDE SLIDER HEARING BY 
C STOCHASTIC APPROACH 
C 

C THIS PROGRAM CONSISTS OF 3 SUBROUTINES 

C 1 - subroutine TABLE IS TO FORM ThE T1* T2» T3 ARRAYS WHICH WILL BE USED TO 

c COMPUTE The pressure profile 

C 2 - subroutine FUNC is to evaluate the functions Fl. F2 AND F3 

C 3 - subroutine INT is ThE INTEGRATION ROUTINE TO CALCULATE THE PRESSURE 

C PROFILE AND then ThE TOTaL LOAD 

C The OUTPUT OF THIS PROGRAM ARE ■ wLOADI » WLOAD2 » WL0AD3 AND THE CORRESPONDING 
C OSIG(I) 

C WLOADI = LOAD IN SLIDING ROUGH SUPFACEtFIXEO SMOOTH SURFACE CASE 

C WL0AD2 = LOAD IN FIXED RCLGH SURF ACE »SL IDING SMOOTH SURFACE CASE 

C WLOAD3 = LOAD IN BOTH SURFACES WITH SAME ROUGHNESS DISTRIBUTION 

C 

C DATA CARO NO. 1 
C KF = total no. of GRID POINTS 
C NF = NO. OF DSIGd) DATA 

C PI = 3.1A1S>»3 

C DATA CARD NO. 2 

C DSIG(I) = DATA FUR SIGMA / H(MIN) RATIO 
C DATA CARD .'.0.3 

C these data are THE OUTPUT FROM PkqGRAM SUCHAS PROGRAM DATP0L2 WHICH EVALUATE 
C THE integral OF ( DISTRIoUTION FUNJCT I 0N*HSS ) / 11 . +SH'->HS5 ) **3 AND 

C (DlSTRiBUTION F UNCT I ON/ ( 1 . ♦ SH*hSS > **3 

C TnIS DISTRIbUTlON FUNCTION CAn BL POLYNOMIAL .GAUSSIAN, SINUSOIDAL OR ANY 
C OTHER FUNCTION 

C DSHfl) = THE ABSCISSA RANGING FROM 0 TO 0.333 
C UV(I) = THE NONDiMENSlONAL DSHII) 

C 064(1) = THE INTtGRAL Of 35 ( 1 .-hSS*»2/9) *»3/ ( 1 -.SH*HSS) **3*HSS 
C IN The case of polynomial distribution 

C U62(I) = THE InTLGRAL OF 35 /96* ( 1 .-HSS«*2/9 J *«3/ ( 1 +SH*HSS ) **3 
C IN THE CASE OF POlYNCMIAL DISTRIBUTION 

C CG5(I) = 0G2(I)/u64(l) 

f 

C 

common DSh(201 ) ,DV (201) .CC-4 (201) . UU5 (20 1 ) . DG2 ( 20 1 ) .WLOADI (14) , 

1 WLOAC«: (14 ) .WLUA03 (14) ♦DSIG(14),Tl(l4.i01).T2(l4,101).T3(l4.101)r 

2 KK.KF.NF.S.x.F 1.F2,F3.C1.C2.CJ.PI.XX(101) 

PRINT 1 

read 2.KF.NF.PI 
print A.KF.NF.RI 
READ 3. (DSIGd). 1 = 1. NF) 

PRINT 4, (DSlG(l) »I=1.NF) 

DO 100 1=1.201 

100 READ b.DSHd ) .DV (1 ) .DG4(1) .UGbd) .062(1) 

DO 99 J=1.KF 
99 XX(J) = (J-l)*0.01 

WREF = ALCG(2.) - 2./3. 
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DO 101 I=1»NF 
KK = I 

S = DSIG(KK) 

CALL TABLE 

Cl = (-TKKKfKF.) + S*T3(Krt»KF)) / T2(KK,KF) 

C2 = <-Tl(KK»KF) -S*T3(KK*KF) ) / T2(KK»KF) 

C3 = -T1(KK»KF) / T2(KK*KF) 

CALL INT 

101 CONTINUE 

PRINT ®,wRer 

OC 102 I = l» Np 

PRINT /»I»DSlG{l) * wlOAU 1 (i) »rtLuAD2(I) .WL0AD3(I) 

102 continue 

c 

1 FORMAT (l Hi /2 X»^»kF*»3X »*nF«» 1 ox* *PI*/) 

2 format (2I5*F20.1A) 

3 format (20F5. 3) 

A FORMAT (X,'-»(DSIG( I) . 1 = 1 *NF ) <»*bX*20 (X.F5.3)/) 

5 FORMAT (F10.6*Fi0.2* JF2O.10) 

6 FORMAT (///X*»wHtF = <^*El7,lC/3X»*I**6X»«uSIG**12X**lKL0ADl*»14X* 

6 »wL0Au2<^* 1AX**wlOAC3«/) 

7 FORMAT ( I5*5x*F3.3*b(3X*E17 .10) ) 

STOP 

E N ij 

SUbROUl INE TABi-t 
C 

COMMON DSH (201 ) *DV (20i ) . CG A ( 201 ) » 0G5 ( 20 1 ) *DG2 (201 ) *WL0AD1 (14) * 

1 WLOAD^: (lA) .WLUAD3 (lA) » OS I G ( 1 A > , T 1 ( 14 » 1 0 1 ) * T2 ( 1 A * 1 0 1 ) * T3 ( 14* 10 1 ) » 

2 Kn,KF*NF*S,X.H*F2,F3*C1*C2*C3,PI,xX(101) 

c 

DIMENSION El (iOl) *E2(i01) *E3(i01) 

J = 1 
X = XX(J) 

CALL FUNC 
El (J) = FI 
E2(J) = F2 
E3(J) = F3 
T1(KK*J) = 0. 

T2(KK,J) = 0. 

T3(KK,J) = 0. 

DO 201 J =2.KF 
X = XX (J) 
call FUNC 
El ( J) = FI 
E2(J) = F2 
E3(J) = F3 

X = (XX ( J-1 ) aXX ( J) )/2. 

CALL FUNC 

T1(KK»J) = T1(NK*J-1) ♦ C.01/2./3.*(El (J-i) ♦ A.«Fi ♦ EKJ)) 

T2(KK,J) = T2(M<*J-1) + 0.01/2./3.*(E2(J-i) + A.*F2 ♦ E2(J)) 

T3(KK.J) = T3(NK*J-1) ♦ 0.01/2./3.*(E3(J-l) ♦ A.*F3 ♦ E3(J)) 

201 CONTINUE 
RETURN 

END 
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SUbROtillNE FUNC 

common DSH(20D »0V <20 1 ) ♦CG** (201) ♦065(201) »DG2(20i ) ♦WLOADl (14) ♦ 

1 WLOAD-t: (14) ♦WL0AD3(14) tOS IG ( 14 ) . T1 1 14 ♦ 1 0 1 ) ♦ T2 ( l4 ♦ 1 01 ) ♦ T3 ( 14» 101 ) t 

2 KK*KFfNF^S^XfFl^F2fF3fCltC2tC3tPI^XX(101) I 

H = 2. - X 

Sh = S/H 

V = 1 ,+SH/O. 001665 
MV s V 
VM X MV 

D2 X (AL0G(DG2(MV*1) )-AUCG (0G2(MV) ) ) * (V-VM) ♦ ALOG (DG2 (MV) ) 

62 X EAP(D2) 

04 X (AL0G(-0G4(MV*1) )-ALCG(-DG4 (MV) ) )*(V-VM) ♦ ALOG(- 0G4(MV)) 

G4 X -cXP(D4) 

FI = GC/H«*2 
F2 X Gii/H**3 
F3 = Gh/H**3 
RETURN 

END 

SUbROUilNt InT " ' 

COMMON DSh (201 ) ♦OV (201 ) *DGa ( 201 ) ♦065(201) ^062 <201 ) ♦WLOAOl (14) ♦ 

1 WLOAC^: (14) ♦WL0AD3( lA) ♦USlo(l‘*),Tl(l4^101)^T2(14^i01)^T3(14^101)^ 

2 KKfKF^NF^S.X.Fl^F2^F3fCl^C2^C3.PI^XX(101) 

OIMENSiON Pi(i‘»^10i)^Pi;(lM^101)^P3(14^101) 

PRINT 31 fKKfOSlG (KK) 

00 301 J=1«KF 

PKKK^J) = Tl(M^^j) ♦ C1«12(XK^J) - S*T3(KK^J) 

P2(KK^J> = TK^K^J) ♦ C2«T2(KK^J) ♦ S*T3(KK^J) 

P3(KK,J) X T1(^K,J) ♦ C3«12(KK^J) 

PRINT j 2^ J^Pl <NK^ J) ^P2 (KK^J) ♦?3(KK^J) 

301 CONTINUE 
WLUADl (KK) X 0. 

WL0A02(KK) X 0. 

VVLOAD3 (KK) = 0. 

KFF X ^F-2 

DC 302 J=1^KFF^2 

wLOADUKK) X wLOADl{K^) ♦ 0.01/3.« (PI (KK^ J) + 4 . *P 1 (KK , J* 1 ) ♦ 

1 (KK« J + 2) ) 

WLOAD21KK) = rtl-0AD2(KK) ♦ 0 . 0 1/3 .* (P2 ( NK ♦ J ) ♦ 4.oP2 (KK^ J^l ) ♦ 

2 (KK ♦ j + 2 ) ) 

y»LUAD3lKK) X Wl0AD3(KK) ♦ 0 . 0 1 /3 . * (P3 (KK ♦ J > + 4.*P3 <KK^ J*1 ) ♦ 

3 P3(KKtJ*2)) 

302 CONTINUE 
PRINT a 

PRINT 7^KK,osiG(KK) ♦iVLOACl (KK) ♦WLCAD2 (KK> ♦WL0AD3(KK) 

C 

31 format (/X^«KK**^ I4,5X^<»USIGx*,F 7.a//3X^*J*^ 10a.,*P1*,18X^*P2*^18Xi 
C *P3*/) 

32 F0RMAT(16^3(3x^tl7,10) ) 

6 FORMAT (/3X *<*I**6X^*CS 16* ♦ 12X^* WLOAOl* ♦ 14X^*WL0 a02*^ 14X^*»(L0A03*/) 

7 format ll5^5x»Fa,3^b(3X^E17.10)) 

C 

RETURN 

CND 
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OOOOvOCiuOOOOv'OOCCuOjvU End of heco^<u 
101 10 3.1^1'3'»i:to53‘j9 

0.06 0.1 0.14 O.IB i).22 0.2fe 0.3 0.31 0.32 0.33 

C data for DSH(I>» UV(I), DGh(I). 4G2(I)» DG5(I» ARE OMITTED 
OOOOOOOOOOOOuOOOOOOOwO END OF INFORMATION 
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PROGRAM DATPOLc ( iNPuT .OUPUT*PUNCH) 

C 

C THIS program is to calculate CaTA For functions G2 and ga for The case when 
c The asperity height list-'Icut ion function is in the form of polynomial 
C function 35/96*a-HSS>»*e:/S >**3 

C G2 = THE INTEGRAL OF Gb/Uftc ( i-nSS**2/9 ) **3/ ( 1 *SH*HSS ) **3 FOR HSS = -3 TO 3 
C G4 = THE INTEGRAL UF 35/9e« < I -nSS**2/9 ) **3/ ( 1 +Sh*HSS ) *«3*HSS FOR HSS = -3 TO 3 
C SINCE HSS RANGES FnCM -3 TO 3 * Sh MUST BE LESS THAN 1/3, THUS SH VALUES 

C USED IN this program RAnuES from 0 TO 0.333 
C SH<I), I = lf201* WlTn 200 lMFORM GRlOS OF SIZE 0,001665 
C V(I) = NONOIMENSIONAL ShO)' 

c 

c 

C THIS program consists OF THREE SUbROUTiNES 

C 1 - subroutine main is to calculate RII and RI2 FROM SUBROUTINES INTI 
C AND INT2 

C 2 - subroutine IMI IS TO INTEGRATE THE FUNCTION 
C ( 1 ,-HSS**2/9. ) **3/ ( 1 . ♦Sl-«hSS) **3*hSS 

C 3 - SUBROUTINE INT2 IS TO INTEGRATE THE FUNCTION 
C ( 1 ,-HSs«*2/9. ) **3/ ( i .♦Sh*(-SS) **3 

C 

C THE RESULTS OF G2n)»GA(I), GB(I) TOGETHER WITH THE CORRESPONDING SMI)tV(l) 
C values ^HICH ARC THE OuTHlT OF THIS PROGRAM ARE IN THE FORM OF DATA CARD 
C 

c 

COMMON A,B,RI1*RI2,SIMPSCN 

dimension G2 (201) ,0^1(201) *65(201) , V (201 ) , SH (202 ) 

PRINT 1 
SHd )=0.0 
DO 10 1=1*201 
CALL main (SH(D) 

G2(I) = 35. /96. *RI2 
6A(I) = 35. /96. *RI1 
65(1) = 62(1) /Gh(I) 

V(I)=1.0+SH(I)/0.001665 
PRINT J 

PRINT I *SH( I ) * V ( 1 ) ,GA ( I) *G5 (1) *G2(I ) 

SH(I + 1)=SH(D*0. 001665 
10 CONTINUE 

PRINT c 
00*20 1 = 1*201 

PRINT s*SH ( I ) * Y ( I ) *GA( I ) *CB (I ) *G2(1) 

PUNCH 5*Sn(I) *V(I) ,G4(1) *05(1) *G2(I) 

20 CONTINUE 

format (IH l ) 

format (5X**SH**8X**V** 16a**GA** 1BX**G5** laX**G2*/) 

FORMAT (/2x**I-tt*5X»«SM(I)«*5X»<»w(l ) **7X**OA ( I ) ** 15 X**g5 ( I )** 15X* 

C*62( !)<*/) 

FORMAT (X* I3*2X*F9.b*2X*Ke. 1*3 (2X*E17.10)/) 
format (FI 0.6*F10.2.3F20. 10) 

STOP 
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SUoKOLilNr: MAlN(iH) 

COMMON A*B»RI1 »RI2*SIMPSCN 
DIMENSION D(3)» L<3)* bM (3 ) » SIM(3) 

data (D(I) *I=i»2)/-3.» (E(I) ♦I=1*2)/-2.6»3./ 

C 

DO 100 J=l*2 

A=0<J) 

B=E ( J) 

CALL INTl(SH) 

SM(J)=bIKPSON 
CALL INT2(SH) 

SIM( J)=SIMPSON 

PRINT il»J»A,B»SM(J) ,SlM(s.) 

100 CONTINUE 

RIl = SM(1) ♦ SM(2) 

PI2 = SIM(I) ♦ SIM(2) 

PRINT i2»RH»RI2 
C 

11 format II ♦ 4X»»A=*,Fb.3» AX ♦*b=*»F6.3«4X»*SM=*»E17.10»4X* 
C«sIM=*»E17.10/) 

12 FORMAT (X»«RI1=**E17.10«5X»*RI2=*,E17. 10/) 

RETURN 

E>D 

bUUROLIINt INTi(bH) 

COMMON A,b,RIl»RI2»SlMPSCN 

PRINT 21 

HSS=A 

FA=( 1 .o-HSS»«2/S.0)-“«3/ (1.0+SH«HSS)*«3*HSS 
hSS=B 

FB=(1 .l--HSS*«2/9.0) ««3/ (1 .0 + Sm^HSS)**3*HSS 

Fl=FA+t- B 

F2=0.0 

HSS=(A*B)/2.0 

F4=< 1 ,u-HSS*»2/9.0)«*3/ (1 .C+SH*HSS)**3*4.0«HSS 

N=2 

T= (B-4) /N 

SIMPS0-x=T/3.0<* (F1*F2*F4) 

DO 10 ^=l♦500 

T= (B-A) /N 

TN=T/2.0 

F2=F2 + t-'»/2.0 

F^=0.0 

00 20 J=1.N 

HSS=A-TN'*J*T 

Fa=(i,Li-HSS**2/9.0)««3/(i.U*SH*HSS)»*3»4.O*HSS ♦ F4 
20 CONTINUE 

SIMP=S1MPS0N 

bIMPS0N=TN/3.0* (FI +F2+FH) 

OELTA = /‘BS( (SIMRSON-SIMP)/ ( blMPbQN + SlMP ) ) 

PRINT 22 »N»SIMP» DELTA 
N=2»N 
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IF(DELTA.LE.0.w005) go TC c 6 
10 CO.-MTINUE 

26 PRINT c2»N*SIMPS0N 

C 

21 FORMAT 14X»*N*»oX»*SIMPS0N'->»10a»«DELTA«/) 

22 FORMAT (X»I4»AX*E17,iO»‘+A»Fa.S) 

RETURN 

END 

SUtr-ROullNE INTtlah) 

COMMON A*B»RI1*RI2*SIMPSCN 
PRINT 21 
HSS = A 

FA=(1 ,v;-HSS«*2/9,0) **3/ ( 1 , 0 + SH*HSS ) **3 
HSS = B 

FB=( 1 .o-HSS*«2/9.0)**3/ ( 1 . 0 +SH*HSS ) **3 

F1=FA+FB 

F2=0.0 

hSb= (A+E) /2,0 

F‘+= ( 1 ,u-HSS<»*2/9.0 > **3/(1 . 0 ♦SH*HSS ) . 0 

N = 2 

T= (B-A) /N 

SIMPSCN = T/3.0« (FI +F2+F4) 

DC 10 r>=l»500 

T= (b-A)/N 

TN=T/2.0 

F2=F2+FA/2.0 

Fa=0.0 

DO 20 J=1»N 

HS5=A-TN+ J*T 

F4= ( 1 .C-HSS*«2/9.0 ) **3/ ( 1 .0 + Si-i<»HSS) *<^3*4.0 + FA 

20 CONTINUE » 

SIMP=S1MPS0N 

SIMPS0N=TN/3. 0* (FI ♦F2*F4) 

DELTA=aBS ( (SIMKSUN-blMP) / ( 5 I MR60N + S I MP ) ) 

PRINT i:2»N»SIMR»DELTA 
N = 2*N 

IF (DEL rA.LE.O.UOOtS) GO TC 26 
10 CONTINUE 

26 PRINT «:2»N»SIMPSUN 

C 

21 FORMAT (AX .«N* » bX ♦ *S I MPbON » 1 OX » «DElTA«/ ) 

22 FORMAT (X» I4»AX»E17. 10«ha»F6.B) 

RETURN 

end 

CCO0O00C000000000000v.a end of record 
V vcO'..0u00000u00</000 j wO EisD Ur INFORMATION 


Data 

for G 2 , 

G^ and G^ with polynomial distribution of 

the asperity 

height. 




SH 

* ct/H 




V 

= non-dimensional SH 




- T- 

_ j.* 



°2 

’L[i 

+ 6 (ct/H)J-" 



p 

- r~ 





“L[i 

- 'i 

+ 6 (ct/H)J 



^5 

= G^/G^ 




* 

g(6 ) 

f96l 

9 / 

-3 s 6* s 3 





elsewhere 


SH 

V 

G4 

G5 

G2 

O.QODOOO 

1.00 

. 000 000 0008*7350 8622.89722U4263 

.9999937040 

.001665 

2.00 

-.0049951360 

-200. 1968247935 

1.0000103701 

• 003330 

3.00 

-.0099912933 

-100.0931790408 

1. QQ006Q3Q55 

.004995 

4.00 

-.0149684332 

-66.7276896773 

1.0001435203 

.006660 

5.00 

-.0199876140 

-50.0439939147 

1.0002600317 

.008325 

6.00 

-.0249895170 

-40.0331812483 

1. 0004098635 

.009990 

7.00 

-.0299948249 

-33.3588560465 

1.0005930462 

.011655 

8.00 

-.0350042216 

-28.5911118930 

1. 0008096174 

.013320 

9.00 

-.0400183928 

-25.0149881471 

1.0010596214 

.014985 

10.00 

-.0450380259 

-22.2332815146 

1.0013431094 

.016650 

11.00 

-.0500638108 

-20.0076686839 

1.0016601395 

.018315 

12.00 

-.0550964396 

-18.1864886970 

1. 0020107769 

.019980 

13.00 

-.0601366075 

-16.6686338844 

1.0023950934 

.021645 

14.00 

-.0651850124 

-15.3841064170 

1. 0028131682 

.023310 

15.00 

-.0702461021 

-14.2821460160 

1. 0032650874 

.024975 

16.00 

-.0753129297 

-13.3277373168 

1. 0037509442 

.026640 

17.00 

-.0803901121 

-12.4924672077 

1.0042708389 

.028305 

18.00 

-.0654783624 

-11.7553127078 

1. 0048248792 

.029970 

19.00 

-.0905783984 

-11.0999222553 

1.0054131799 

.031635 

20.00 

-.0956909428 

-10.5133864718 

1. 0060358632 

.033300 

21.00 

-.1008167234 

-9.9853776741 

1.0066930586 

.034965 

22.00 

-.1059564734 

-9.5075352232 

1. 0073849032 

.036630 

23.00 

-.1111109320 

-9.0730184962 

1. 0081115416 

.038295 

24.00 

-.1162808444 

-8.6761764664 

1. 0QB8731260 

.039960 

25.00 

-.1214669623 

-8.3122998848 

1.0096698164 

.041625 

26.00 

-.1266700439 

-7.9774329347 

1.0105017804 

.043290 

27.00 

-.1318908551 

-7.6682283475 

1. 0113691938 

.044955 

28.00 

-.1371301688 

-7.3818347149 

1. 0122722402 

.046620 

29.00 

-.1423887658 

-7.1158079458 

1. 0132111113 

•048285 

30.00 

-.1476674354 

-6.8680410430 

1.0141860071 
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.049950 

31.00 

-.1529669751 

- 6.6367079229 

1,0151971358 

.051615 

32. 00 

-.1582881916 

- 6.4202181087 

1. 0162447141 

.053280 

33.00 

-.1636319007 

- 6.2171799186 

1, 0173289673 

.054945 

34.00 

-.1689989282 

- 6.0263703463 

1. 0184501295 

.056610 

35.00 

-.1743901098 

- 5.8467102566 

1.0196084435 

.058275 

36.00 

-.1798062917 

- 5.6772438332 

1. 0208041609 

.059940 

37.00 

-.1852483314 

- 5.5171214518 

1.0220375429 

.061605 

38.00 

-.1907170973 

- 5.3655853297 

1. 0233088596 

.063270 

39.00 

-.1962134701 

- 5.2215363317 

1. 0245357631 

.064935 

40.00 

-.2017383425 

- 5.0802341398 

1. 0248780151 

.066600 

41.00 

-.2072926201 

- 4.9508585771 

1.0262764462 

.068265 

42.00 

-.2128772216 

- 4.8277327210 

1. 0277143282 

.069930 

43.00 

-.2184930795 

- 4.7104100548 

1. 0291919985 

.071595 

44.00 

-.2241411406 

- 4.6027585049 

1.0316675410 

.073260 

45.00 

-.2298223663 

- 4.4957019781 

1. 0332128666 

i 074925 

46.00 

-.2355377334 

- 4.3933455139 

1.0347986442 

.076590 

47,00 

-.2412882345 

- 4.2953824033 

1. 0364252365 

.078255 

48.00 

-.2470748786 

- 4.2015320339 

1.0380930172 

,079920 

49 . 00 

-.2528986917 

- 4.1115371720 

1. 0398023717 

.081585 

50.00 

-.2587607174 

- 4.0251615764 

1.0415536970 

.033250 

51.00 

-.2646620173 

- 3 . 942 ia 7 RP 97 

1.0433474022 

.034915 

52.00 

-. 27 Q 6036722 

- 3.8624156356 

1. 0451839087 

,096580 

53.00 

-.2765867821 

- 3.7856604805 

1. 0470636503 

.038245 

54.00 

-.2826124672 

- 3.7117508799 

1.0489870739 

,089910 

55.00 

-.2886818687 

- 3,6405287375 

1.0509546391 

.091575 

56.00 

-.2947961494 

- 3,5718472638 

1. 0529668195 

,093240 

57.00 

-.3009564942 

- 3.5055701487 

1. 0550241022 

.094905 

58.00 

-.3071641114 

3.4415706427 

1. 0571269883 

.096570 

59.00 

-.3134202331 

- 3.3797307319 

1. 0592759937 

,098235 

60.00 

-.3197261160 

- 3.3199403982 

1. 0614716488 

. 0-99900 

61.00 

-.3260830426 

- 3,2620969524 

1. 0637144994 

,101565 

62.00 

-.3324923217 

- 3.2061044340 

1. 0660051069 

.103230 

63.00 

-.3389552897 

- 3.1518730679 

1. 0663440488 

,104895 

64.00 

-.3454733110 

- 3.0993187738 

1. 0707319187 

.106560 

65.00 

-.3520477796 

- 3.0483627208 

1. 0731693274 

.108225 

66.00 

-.3586801197 

- 2,9989309240 

1. 0756569029 

.109890 

67.00 

-.3653717869 

- 2.9509538769 

1.0781952911 

.111555 

68.00 

-.3721242692 

- 2.9043662172 

1.0807851561 

.113220 

69 . 00 

-.3789390883 

- 2.8591064219 

1. 0834271808 

.114885 

70.00 

-.3858178006 

- 2.8151165293 

1.0861220677 

.116550 

71.00 

-.3927619984 

- 2.7723418840 

1, 0888705387 

.118215 

72.00 

-.3997733116 

- 2.7307309043 

1. 0916733368 

.119880 

73.00 

-.4068534083 

- 2.6902348687 

1,0945312255 

.121545 

74,00 

-.4140039967 

- 2.6508077195 

1. 0974449905 

.123210 

75.00 

-.4212268264 

- 2.6124058831 

1,1004154394 

,124875 

76.00 

-.4285236897 

- 2.5749881038 

1 . 1034434032 

.126540 

77.00 

-.4358964234 

- 2.5385152918 

1. 1065297365 

,128205 

76.00 

-.4433469102 

- 2.5029503817 

1.1096753182 

.129870 

79.00 

-.4508770805 

- 2.4682582034 

1,1128810527 

.131535 

80.00 

-.4584889140 

- 2,4344053612 

1,1161478703 
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.133200 

.134865 

.136530 

.138195 

.139660 

.141525 

.143190 

.144855 

.146520 

.148185 

.149850 

.151515 

.153180 

.154845 

.156510 

.158175 

.159840 

.161505 

.163170 

.164835 

.166500 

.168165 

.169830 

.171495 

.173160 

.174625 

.176490 

.178155 

.179820 

.181485 

.183150 

.184815 

.166480 

.188145 

.189810 

.191475 

.193140 

.194805 

.196470 

.198135 

.199800 

.201465 

.203130 

.204795 

.206460 

.208125 

.209790 

.211455 

.213120 

.214785 


81.00 

62.00 

83.00 

84.00 

85.00 
86 . 00 

87.00 

88.00 
89 . 00 

90.00 

91.00 

92.00 

93.00 

94.00 

95.00 

96.00 

97.00 

98.00 

99.00 
100.00 
101.00 
102.00 

103.00 

104.00 

105.00 

106.00 

107. 00 

108.00 

109.00 

110.00 
111.00 
112.00 

113.00 

114.00 

115.00 

116.00 

117.00 

118. 00 

119.00 

120.00 
121.00 
122.00 

123. 00 

124. 00 

125.00 

126. 00 

127.00 

128.00 

129.00 

130.00 


-.4661844416 
-.4739657472 
-.4818349697 
-.4897943053 
-.4978460090 
-.5059923974 
-.5142358509 
-.5225788159 
-.5310238076 
-.5395734124 
-.5482302910 
-.5569971807 
-.5658768992 
-.5748723473 
-.5839865120 
-.5932224706 
-.6025833939 
-.6120725501 
-.6216933087 
-.6314491450 
-.6413436442 
-.6513805062 
-.6615635504 
-.6718967209 
-.6823840917 
-.6930298724 
-.7038384143 
-.7148142166 
-.7259619328 
-.7372863781 
-.7487925362 
-.7604855676 
-.7723708173 
-.7844530238 
-.7967403261 
-. 8 Q 92362831 
-.8219478646 
-.8348814812 
-.8480437364 
-.8614416905 
-.8750623736 
-.8889732995 
-.9031222298 
-.9175372398 
-.9322267347 
-.9471994673 
-.9624645568 
-.9780315083 
-.9939102345 
- 1.0100755107 


- 2.4013601233 

- 2.3690923186 

- 2.3375732418 

- 2.3067755645 

- 2.2766732533 

- 2.2472414935 

- 2.2184566178 

- 2.1902960406 

- 2.1627381961 

- 2.1357624808 

- 2.1093491996 

- 2.0834795170 

- 2.0581354079 

- 2.0332996159 

- 2.0089556112 

- 1.9850875521 

- 1.9616802495 

- 1,9387191326 

- 1.9162541838 

- 1.8941383961 

- 1.8724283612 

- 1.8511116757 

- 1.3301764135 

- 1.8096111019 

- 1.7394047000 

- 1.7695465774 

- 1,7500264946 

- 1.7306345845 

- 1.7119613347 

- 1.6933975710 

- 1.6751344411 

- 1,6571634006 

- 1.6394761979 

- 1.6220648613 

- 1.6049216864 

- 1,5880392234 

- 1,5714102663 

- 1.5550278414 

- 1.5388851972 

- 1.5229757946 

- 1 . 507293297 C 

- 1.4918315618 

- 1,4765846316 

- 1.4615467257 

- 1.4467122327 

- 1,4320757029 

- 1.4176318408 

- 1.4033754984 

- 1.3893016685 

- 1.3754539099 


1,1194767281 
1. 1228686109 
1.1263245322 
1. 1298455350 
1. 1334326929 
1. 1370871108 
1, 1408099265 
1, 1446023114 
1. 1484654718 
1,1524006499 
1.1564091255 
1, 1604922170 
1. 1646512829 
1.1688877229 
1. 1732029801 
1. 1775985421 
1. 1820759426 
1, 1866367633 
1, 1913224038 
1. 1960520706 
1.2008700286 
1.2057780603 
1, 2107780060 
1.2158717656 
1, 2210613009 
1, 2263486387 
1.2317358729 
1,2372251675 
1.2428187595 
1,2485189617 

1.2543281667 
1.2602488493 
1,2662835709 
1.2724349830 
1.2707058309 
1.2850989587 
1.2916173126 
1.2982639475 
1, 3050420294 
1, 3119548430 
1.3190057961 
1, 3261984258 
1. 3335364049 
1.3410235485 
1,3406638207 
1.3564613430 
1.3644204014 
1. 3725454555 
1. 3808411471 
1.3893123104 
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.216450 

131.00 

- 1.0266066512 

- 1.3617328315 

1. 3979639819 

•218115 

132.00 

- 1.0434818142 

- 1.3481800950 

1.4068014115 

.219780 

133.00 

- 1,0607127742 

- 1.3347911976 

1.4158300742 

•221445 

134.00 

- 1.0783118426 

- 1.3215617472 

1.4250556828 

.223110 

135. 00 

- 1.0962918997 

- 1.3084874581 

1. 4344842012 

.224775 

136.00 

- 1.1146664290 

- 1.2955641449 

1.4441218590 

•226440 

137.00 

- 1.1334495539 

- 1.2827877178 

1.4539751665 

•226105 

138. 00 

- 1.1526560771 

- 1.2701541777 

1. 4640509318 

.229770 

139.00 

- 1.1723015235 

- 1.2576596113 
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.299700 181,00 -2. 8<»4857<i004 
.301365 182,00 -2.9295837773 
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Data for G2 and with sinusoidal distribution of 
SH *= a/H 

V non-dimensional SH 
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,91499997 
,23177560 
,58087237 
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PROGRAM RCLWAVI ( INPUT *OL IPUT ) 


C 

C THIS program is to COHPuTe Trt6 EFFECT OF SURFACE ROUGHNESS AND WAt/INESS ON THE 
C iNTEGRATtD PRESSURE OF AN EhO CONTACT 
CONLY THE CASE OF pU^E ROLLING IS STUDIED 

C ASSUMING SINUSOIDAL wAvlNESS AND SINUSOIDAL ROUGHNESS DISTRIBUTION 
C 

C THIS program CONSISTS OF 3 SUBROUTINES 

C 1 - s'ubROUUNE GRUtsIN IS TC CALCULATE WGRUB • THE REDUCED PRESSURE AT THE 
C INlET by the smooth film GRLBIN APPROACH 

C 2 - SUBROLflNE ROUGH IS TO CALCULATE WROUGH(I) t THE REDUCED PRESSURE AT 

c The Inlet by stochastic ihecky for rough surfaces 

C 3 - subroutine WAVI IS TO CALCULATE WSlN<I) » THE REDUCED PRESSURE AT THE 
C INLET FOK THE SINUSOIDAL wAViNESS SURFACE PROFILE 
C 

C DATA CARD NO. 1 

C KF = TOTmL number OF GRID POINTS 

C KM = GRiU NUMBER * XX (KM) = -2,0 

c NAMP = nomber of AHFLITUUE CATA 

C NCYC = NUMBER OF WAVINESS CYCLES DATA 

C data CARD NO. 2 
C HO = center FILM THICKNESS 
C w = OimEhSiONLESS CONTACT LCAD 
C PI = 3.1‘»1593 
C DATA CARO NO. 3 

C OCYC(I) = waviness CYCLES data 
C DATA CARD NO. 4 
C DAMP(l) = AMPLITUDE DATA 
C DATA CARO NO. 5 

c These data are the output from program suchas program datpol 2 which evaluate 
C the IMEGRAL OT < distribution FUNCTI0N<>MSS)/a.+SH<‘HSS)*«3 AND 
C (CISTwiBUTION rUNCT 10N/a.+SHttHSS)*«3 

c iHis distribution Function can be polynomial ,gaussian» sinusoidal or any 
C other t unction 

C OSH(I) s the abscissa ranging from 0 TO 0.333 
c UV(I) = THE NONDiMEMSIONAL DSH ( I ) 

C OGh(I) = The INTlGRAL of 3S/96*<l.-MSS**2/9)*»3/(i+SH*hSS)**3*HSS 
C IN THE CASE OF POLYNOMIAL DISTRIBUTION 

C 062(1) = THE INTlGRAL OF 3£ /96* ( 1 . -hSS**2/9 > **3/ ( 1 +SH*HSS ) **3 
C IN THE CASE OF POLYNOMIAL DISTRIBUTION 

C UGSd) = DG2U)/0G4(1) 

C 

C 

COMMON KF»KM,N*MP»NCYC»KPK *kmk»kff 

COMMON HOtWtPI •wGRUB.CYCfAHp.KK.S » DG2 (20 1 ) *064 ( 20 1 ) 
common XX ( 1501 ) »H( 1501) , W£ IN ( lb) , WROUGH ( 1 6) . 6 (ISO 1 ) 

DIMENSION DCYC(20) .DAMP (16) .SwiSINdb) .SWROdb) .PERR(lb) ,OSH(201) , 

! DV(20i) 

C 

print i 

READ 2.KF.KM.NAMP.NCYC 
READ 3, HO, w, PI 
READ 3, (DCYC(i) ,I=1 .NCyC) 
read 3, (DAMP( I ) , I=1,NAMF) 

DO IOC 1=1.201 

100 read 5.0SHd),Dv(I).DG2(I),D6M(I) 
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PRINT si 

PRINT if » KF ♦ KM » KAMP ♦ NCYC 
PRINT b2 
PRINT ‘♦♦H0»W*PI 
PRINT 33 

PRINT '♦t(DCYC(I)*I = l»NCYO 
PRINT 

PRINT (L)AMP(i) »I = 1»NAMF) 

KKM = KM ♦ 1 
KMK = KM - 2 
KFF = KF - 2 
^<X(1) = -5. 

00 101 J =2*KM 

101 XX(J) = XX(J-l) ♦ 0.003 
DO 102 J =KMMfKF 

102 XX(J) = XX(J-l) ♦ 0.002 
DO 103 J=1»KF 

X = XX(J) 

H(J) = HO ♦ 4,<>W/PIi>(AbS(X)*SQKT{AbS( (-X)^»«2-l.) ) - AL06(ABS(X) 
H SORT(mBS( (-X)<Hf2-l. )))) 

103 CONTINUE 
PRINT bS 

PRINT H.XX (KM) »XX (KF) 

C 

call GnUBIN 
DO 104 I=1»NAMP 
KK = I 

S = DAfiP(I) / 3, * HO 
CALL ROUGH 

S\>.RO(I) = WROUGh(I) / ir.GRLB 

104 COi'.TINUE 

DC lOOC 1I=1»NCYC 
CVC = i-CYC(II) 

PRINT 7, IlfCYC 
DO 105 I = I.NAMP 
KK = I 

AMP = uAMP(I) * nO 
CALL l««VI 

SWSIN(l) = WSINd) / wGkLc 
PERP( l) = (WSlN(I) - wkOlGH(I)) / 

PRINT o, l»AMP »WS1N< I ) *SwiiN( I ) »ta6GuO*it I ) ) tPERPn) 

105 CONTINUE 

DC lOfc I=1.NAMP 
KK = 1 

AMP = -DAMP(I) * HO 

call I«aVI 

SwSIN(i) = WSINd) / WbRLe 

PERRd) = (HSINd) - wHOLGhd)) / uROUGhID 

PRINT o. I»AMP »wSiN d ) »SwSlN ( 1 ) ♦t.HOuGt't I » ♦ SwRC ( I ) tPERR d) 

106 CONTINUE 

1000 continue 
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c 

1 FORMAT(lHl) 

2 FORMAT (2015) 

3 FORMAT (8F10. 7) 

A F0RMAT(9(2X»E13,6) ) 

5 FORMAT (2Fi0.6»3F20. 10) 

6 FORMAT (I5»8(2X»E13, 6) ) 

7 FORMAT (////////// 3 X»*I I=«» I 3 t 10x«»CYC=*«F(>. !//<•* «*I«»6X»*DAMP»* 11 X 
7 »*WS1W*» 10X»-»SmSlN*»yX»«»hOU&>^** 1 lX»*i.»Ru*« 11X«*RERR*/) 

51 format (/3X**KF<^»3X»«i<M*»X»*NAMH*»X»*NCyC*) 

52 FORMAT (/6X»*H0*»14X»*ft'*f 12X,*r?l*) 

53 F0KMAT(/6X**DCVC(I)*) 

5a F0PMAT(/6X»«DAMP(I)«) 

55 format (/6X««XX (KM) » »9X»*XX <KF) *) 

c 

STOP 

— 

SUbROUlINt GRUblh 

C 

COMMON KF*KM»NAMP,.NCYC*KMM «KMKtKFF 

COMMON H0*W*PI »WoRUB*C'IC»AmP*KK»S »DG2 (201 ) »0GA (201 ) 

COMMON XX(1501)*H(1501)fWSlN(lt>) ,h ROUGH {16)»G(150i) 

C 

DO AOl J=1»KF 

G(J) = (H(J> - HO) / h(J)«*3 

AOl continue 

wGRUB = 0, 

DC A02 J=1»KMK»2 

A02 WC-RUB = WGRUB ♦ 0.003/3.* (G(J) •» ^.«G(J^1) ♦ G«J*2)> 

DO A03 J=KM,kFF»2 

A03 wGRuB = WGRUB ♦ 0 . 002/3 .«( C- ( J) ♦ A.*G(J*D ♦ G(J*2)) 

PRINT Alf HOfiVGRUB 

A1 FORMAT {/////5X «*HO=*»FiO .c »20X»*WGRUb=**ElG.6//) 

C 

RETURN 

END 

SUbROu I INE RUU'Jh ' 

Common Kr»KM,N'bMp.NCYC.KMM * kmk»kff 

COMMON hO.W.Pl .WG«U6»CYC* AMPtKK.S . DG2 (20 1 > » DGA (20 1 ) 
common XX(i501»*n(1501)*wilN(it>) .WR0UGH(16) *6(1501) 
dimension TR(1j01) 

DC 301 J=1*KF 

V = 1. + 5 / H(J) / 0.001C6S 
MV = V 
VM = MV 

02 = (AL0G(DG2(MV + i) ) - AL OG ( DG2 ( MV ) ) ) * (V-VM) ♦ ALOG (DG2 (MV ) ) 

G2 = £XP(D2) 

301 TP(J) = 6(J) * G2 
WROUGH(KK) = 0-. 

DO 302 J= 1»KMK,2 

302 WROUGI-(KK) = WkOUGH(KK) ♦ 0 ,003/3.* (TR (J) ♦ A.*TR(j4l) + TR(J*2)) 
DO 303 J =KM.KFF,2 

303 WROUGh(KK) = WkOUGH(KK) ♦ 0 ,0u2/3 . * ( TR ( J ) ♦ A.*TR(J*1) ♦ TR(J+2)) 
RETURN 

END 
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SOb^OUriNE WAVi 

c 

COMMON KF*KM»NAMP,NCYC»KKN *KMK»kFF 

Common ho*iv*pi »wgru6»cyc»amp*kk»s »OG2(20i) * dg^(201) 

COMMON XXdSOl) »H(1501) fhSlNCib) »KPOuGH(i(S) »G(1501) 
dimension, T(15U1) 

C 

DO 201 J = 1»kF 

del = Amp * COb(e:.*CYC*PI«XX(J) ) 

T(J) = (H(J) - HO) / (hOv) ♦ DEL)**3 ■ 

201 CONTINUE 
l«(SIN(KN) s 0. 

DO 202 Jsl»KMN*2 

202 WSIN(KN) = WSIN(KK) ♦ 0 . 0 0 3/3 , * ( T ( J) * A,*T(J*1) ♦ ,T(J*2)) 
DO 203 J=KM.KFF*2 

■203 WSIN(KN) = WSIiN(KK) ♦ 0 . 0C2/3. » <T ( J> ♦ 4.*T(J + i) ♦ T<J*2)I 
C 

RETURN 

E^^D 

OOOOuOOOOOOuijOOOOOOOuO End of htCORU 
1501 1001 9 14 

0.00001 0.00003 3.14159265 


1. 

2 • 

3. 

A. 

5. 

10. 

15. 

20. 

25. 

o.is 

30. 

0.3 

35. 

0.45 

40. 

0.6 

45. 

0.75 

50. 

0.9 

0.93 

0.96 

0.99 

OATfi 

FOR 

DSH(I) » DV(I) , 

DG4 ( I ) , 

DG2(1)« DG5<I) 

UJ 

OMITTED 



OOOOOOOOOOOOuOOOOOOOuO END OF INFORMATION 
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PROGRAM ASPERTY ( INPUT *Ol.TPuT ♦PUNCH) 

C THIS PROGRAM IS TO SOLVE FOP T*it PtiPTUPBEO PRESSURE DUE TO A SINGLE 3-0 
C RIGID asperity wIThIN AN EL A?T0HTDPCUYNAMIC CONTACT 
C THIS PROGRAM CONSISTS Of P SU6R0UTINES 

C 1 - subroutine TRAnSFN is LSED TO TRANSORM THE INPUT DATA OF CHENG PROGRAM 
C PROM GLOBAL COORDINATE TO ASPERITY COORDINATE 

C 2 - subroutine INTtRPN IS USED TU INTERPOLATE SOME NEW DATA IN ASPERITY 
C coordinate by SECOND ORDLR INTERPOLATION ROUTINE 

C 3 - SUBROUTINE ABCR IS TO FCkM THE A (2^ ♦ 1 A ♦ ) »B (i;9 ♦ ) ♦ C(29»1A»IA) AND 

C R(29flA) MATRICES 

C A - SUBROUTINE TEFPHI IS TC FORM MATRICES T (29 ♦ 1 A ♦ I** ) ♦ E (30 ♦ 1 A ♦ 1 A ) , F (30» lA) 
C AND SOLVt FOP THL PEHTuRBEC PRESSURE PHI(29^1A) BY THE COLUMNWISE MATRIX 
C INVERSION METHOD 
C 

c 

C DATA CARD NO. 1 

C NCI = NUMBER OF THE MAXIMUM ASPEkITY AMPLITUDE 

C NUB = number of OIFFEREta VELOCITIES DATA 

C NGAMAR = NUMBER OF THE ELLIPSTICITY RATIO OF ASPERITY 
C NPLOT = 0 vOO NOT PLOT ThE RESULTS 

C NPRINT2 = 0 ♦ DO NOT PRINT ThE keSULTS 

C imprint = 0 ♦ DO ImOT print TrL RESULTS 
C 

C DATA CARD ImO. 2 

C BSTAR = THE RATIO OF THE MINOR AXIS OF THE ASPERITY TO THE HALF HERTZIAN 
C CONTACl width 

C X3 = THE DISTANCE BETwEEN THE ASPERITY CENTER AMD ThE CONTACT CENTER 

C DI = THE DISTANCt IN X-OIRECTION WHERE THE PERTURBED PRESSURE PHI IS 

C ASSUMED TC EE AERO 

C 

C DATA CARD mO. 3 

C GAMMA(I) = ELLIPSTICITY kATIO DATA 
C DATA CARD MO. 4 

C]H(i) = asfarity amplitude data 

C DATA CARD NO. 5 
C Ul(I),u2(I) = VELOCITY DATA 

C DATA CARD NO. 6 
C PHZ = HEhTIZIAM hRESSuRE 

C AL3E = the product OF PRtSSLRE VISCOSITY COEFFICIENT AND THE EQUIVALENT 
C YOUNGS MODULUS, 

C data CARO NO. 7 

C HO = CENTER FILM THICKNESS 

C 

C data CARD NO. 8 

C KO = IN oLOBAL COORDINATESf TriE GRID NO. OF HERTZIAN CONTACT CENTER WITH 
C XS(KO) = 0 

C KF = IN GLOBAL CDOkDINATESv ThE ORID NO. OF CONTACT EXITWITH XS(KF) = 1 

C NC = IN »*SPERITY COORDINATE^ THE GRID NO. OF HERTZIAN CONTACT CENTER 

C NXvNY = THE NO. OF TOTAL GRID POINTS IN X AND Y DIRECTION OH ASPERITY COORD. 

C IDEG=2 MEANS A PARABOLIC IN TEhPOL aT ION ROUTINE USED IN SUBROUTINE INTERFN 

C 
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C DATA CARD ixO. 9 

C DXS(I) = t-RlD SUES DATA IN OLOSAU COORD. (FKOH ChENG PROGRAM) 

C data card 'nO. 10 

C PS(I) = STEADY state PRESSURE DATA FROM CHEN& PROGRAM 

c data card (< 0 . 11 

C ►^SU) = STEADY STaTE FILM TFICKNESS DATA FROM Ch£NG PROGRAM 
C DATA CARD i»0. 12 

C OX (I) = ORID SIZES DATA IN X-DIRECTION OF ASPERITY COORO. 

C DATA CARD hO. 13 

C DY(I) = ORIC SIZES DATA IN Y-DlRtCTION OF ASPERITY COORD. 

C DATA CAOD NO. 14 

C MIN(I) = NO. USED IN SUBROUTINE INTERPN TO IDENTIFY THE GRID POSITION 
C 

c 

COMMON BSTAR*X3»C1 
COMMON alpha. ALAmDA* GIN V 

common NX.NXI .NY.NYl .DI.CX (28) * DY ( 14) .DEL (29. 15) . HTEP(29«1S)« 

C H1HTEP(29.15) * Pl(29). Pl(29). X(29) 

COMMON AC29.14.14). B (<:■:«. lA. iS) . C(29.l4.14). R(29»14) 

COMMON PHI (29. lA) .PiPhI (25. lA) *HT (29.1b) . NPR INT . NPR INT2 
C 

dimension DXS(e0).X5(80)./XS(6u).PS(80).Hb(80).MlN(29). 

1XCC(29) .PPL0T(t:9) . Y ( 15 ) . h 1 EP (29 ) . C I H ( 5 ) . EP 1 (29 ) . U1 (3 ) »U2 ( 3) 
dimension GAHA«(5) 

READ constant 

read 3.NC1.NUd.nGAmAR»NFLOT .NPR In T2 .NPR INT 
READ ^.BSTAR.a3.DI 
read c. (OAMAR ( I ) . I=1.NGAMAR) 
read (C1H(1) . I = l.NCi) 

DC 100 1=1. NUB 
100 read j.UKI). U2(1) 

READ £:.PhZ.AL3E 
READ 2*HO 

READ S.KC.KF»NC.NX.NY. ICEG 
KFF = r.F-1 

NX I = lxX-1 
NYi = ixY-l 
NCC = NC +1 
ALPHA = AL3E«PH2 

C read smooth FILM SOLUTION DATA 

READ 2. (DXS(I) .I=1.KFF) 

READ J. (PS( I ) » I=l.KO) 
read 3. (hS(I) ♦1=1.aO) 

READ i:. (DX(I) .1 = 1,NX1) 

READ 2. (DY(I) .I=1.NY1) 

READ =>. (MIN(D .1 = 1.NC) 

PRINT 1 
PRINT S3 

PRINT o.NCl .NUe .NGAM,Ak.NFLDT.NPRINT2»NPRINT 
PRINT 31 

PRINT A.BSTAR.A3.DI 
print *.4 
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PRINT (GAMAR (I ) »I = l»NG/if'AR) 

print 

PRINT *♦* (CIH(I) ♦I = I»NC1) 

PRINT 32 
00 99 1=1 f NUB 

99 PRINT '»»I»U1 (I) ♦U2U) 

print 33 

PRINT ‘‘♦*PHZ»AL3E»ALPHA 
PRINT JA 
PRINT A, HO 
PRINT 36 

PRINT b,KC»KF<NC*NX»NY« ICEG*KFF*NCC*NX1*NY1 

AS(1) = “5. 

00 101 I=2»KF 

101 XS(I) = XS(I-l) ♦ DXS(I-I) 

DO 102 I=1*K0 

102 XXS(I) = {XS(I)-X3)/BSTAR 
X(l) = DI 

DO 103 I=2*NX 

103 X(I) = X(I-1> ♦ OX(I-l) 

DO 106 I=1*NX 

106 XOC<I) = X(I)»bSTAR ♦ X3 
Y(l) = 0. 

DO 110 J=2»NY 

110 Y(J) = Y(J-l) + DY(J-l) 

print 38 

PRINT (I«PS(1) ♦I=ltKO) 

PRINT 39 

PRINT 7* (I»HS(I) ♦!=! ♦KO) 

PRINT aO 

PRINT 7, (I,OXS(I> »I=1»KFF) 

PRINT 37 

PRINT 7, (I ♦XSl I ) ♦I=1»K0) 

PRINT a5 

PRINT 7, (I,XXSd) »I = 1»K0) 

PRINT a2 

PRINT 7, (I»DX(1) »I=1*NX1) 

PRINT A3 

PRINT 7, ( I ,DY ( i ) *I=1»NY1) 

PRINT *+7 

print 7. (IiXd) »I = 1*NX) 

PRINT 30 

PRINT 7, (I*XOCd> d = liNX) 

PRINT 31 

PRINT d (J»Y(J) *J=1»NY) 

PRINT Al 

PRINT 6, (MIN(I) »I = 1»NC) 

C 

DO 10*4 1 = 1 »NC 

lOA CALL THANSFN(XxS»PS*HS»KC»X ( I ) fPl ( I ) ♦Hi ( I ) tMIN( I) » IDE6) 
IF (NC .EG. NX) GO TO 107 
DO 105 1=NCC^NX 

Pld) = SQRT(AbS(l. - X0C(I)**2>) 

105 HI d) = 1. 

107 CONTINUE 
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111 


DO 111 I=1»NX 
EPl(I) = EXP(-P1 (I)*ALPHA) 

HIEP(I) = HI (I)*<^3«EPl U) 

PRINT 48 

PRINT 7* (I*P1 (I) »I=1»NX) 

PRINT 49 

PRINT 7,(I,h1(I)»I=1»NX) 

PRINT 54 

PRINT 7» (i,EPl (I) *I=1*NX) 

PRINT 55 

PRINT 7* (I,h1EP(I) *I=1*NX) 

G 

DO 1000 IIK=1»NGAMAR 
GAMMA = GAMAR(IIK) 

GINV = l./GAMMM**2 
DO lOOu IIJ=i»NUB 
UD = (01 (IIJ)-U2(IIJ) )*0.5 
US = (U1(1IJ)+U2(IIJ))*0.5 
ALAMDA = ‘♦B.*UU/h0**2*BbTA« 

ALAM2 = 48,*US/H0**2*BSTAR 
ALAM3 = 4b,/32./ALAM2 
ALAM4 = ALAM3 / GAMMA 
C 

DO lOOw III=I*NC1 
Cl = ClH(III) 

PRINT 56 

PRINT VtllJ* UU» US* ALAMCA, ALAM2* ALAM3* ALAM4 
PRINT 52 

PRINT "s'*!!!* Cl* PH2* UD* AL3E* BSTAR* GAMMA* X3 

TO transform known data tc new coordinates 

TO INTERPOLATE NEW RESULTS FRON KNOWN DATA AT NEW COORDINATES 
CALL INTERPN ( Y*HlEP*EPl ) 

TO FORM MATRICES A(I*J*K)* B(I*J*K)* C(I*J*K)* R(I*J) 

CALL ABCR 

TO FORM MATRICES T(I*J*K)* E(1*J*K)» F(I*J) AND 
SOLVE (■ OR PHId*J)* PlRt-I(I*J) 

CALL TEFPHI 
C 

J = 1 
PRINT B*J 
DO 108 I=1*NX 

108 PPLOT(l) = PHI (I*J) 

C STPLTi 15 A LIBRARY SUBROUTINE TO PLOT RESULTS 
CALL STPLTI ( 1 * A * PPLOT * 29 * 1 M* * 5 * 5HPhI-X ) 

CALL STPLTI ( 1 * aCC»PPLOT*29 * I H* * 7 * 7hPHI-X0C ) 

DO 109 I=i*NX 

109 PPLOT(I) = PlptiI(I*J) 

CALL STPLTI (l*X*PPLCT*29*ih<t,7*s7HPlPHl-X) 

CALL STPLTi (1*XOC*PPLOT*29*1H**9*9HPiPMI-XOC) 

1000 CONTINUE I 
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c 

1 FORMAT tlHl) 

2 FORMAT (8F10. 6) 

3 F0RMAT(5(2X,E13,6) ) 

4 FORMAT (9{2X.e13. 6) ) 

5 FCRMAT(16I5) 

6 FORMAT (18(2X.Ib) ) 

7 FORMAT (8(X»I2»X,E13. 6) ) 

8 FORMAT (////6X»^»J=** 13////) 

9 FORMAT ( I5t8(2X»E13. 6) ) 

10 format (10X»2F10, 6) 

11 FORMAT (X»14E9.<f:) 

31 FORMAT (7X»*BSTAR<^,13X»<‘aJ« » 13X»*D1*) 

32 FORMAT 14X»*I*. /X»*Ui ( 1 ) » 1 IX » *U2 ( 1 ) * ) 

33 FORMAT (7X**PhZ*» 1 1X»<>AL3E« » 1 OX »* ALPHA* ) 

34 FORMAT (//7X»*H0*/) 

35 FORMAT (6X**01HlI)*) 

36 FORMAT (4X»*K0*»bX»«KF**bX »«nC*»5X»«NX*»SX»*NY*»3x»«I0EG*»4X**kFF«» 
6 4A«*^CC*»4X.*NX1«■^^X^^^NY 1*) 

37 F0RMAT(//6X.*X5(I)*/) 

38 FORMAT l//6X,*Pb(l)*/) 

39 F0RMAT(//6X.*hS(I)*/) 

40 format <// bX »*DAS ( 1) */) 

41 FORMAT <6X»*MIN<I)<*) 

42 F0RMATI//6X,*Da(I)*/) 

43 format (// 6x,*Dy(I)*/) 

44 FORMAT (6X»*C'AMMA ( I ) «) 

45 F0RMAT(//6X,*XaS(I)*/) 

47 FORMAT (//6X»*X(I)*/) 

46 Format (// 6x**pi (I)*/) 

49 FORMAT (//6X.*hl (I ) */) 

50 FORMAT 1//6X^*aUC ( I ) */) 

51 format (//6A^*Y (J) */) 

52 format (//2X,*IU*»7X»*C1«» 12X » *PHZ* » 14X » *UD* » 14X » *G* ♦ 1 OX ♦ *8ST AR 

2 *♦ lOX **GAMMA*» 12X»*X3*/) 

53 FORMAT (4X » *NC1 * » hX » *NU b* * X * «NGAm A«* » 2X » *NPLOT* » X * *nPR INT2**X ♦ 

3 *(vPRIi\T«) 

54 FORMAT l//6X»*EHl ( I )*/) 

55 FORMAT (//6X.*hl£R(I)*/) 

56 format (////////// 2X,»1 IJ*»dX»*UD*f 13X»*US**9X»*ALAM0A*t 10X«*ALAM2* 
6 » 10 X*^‘AlAM3»* 10X»*ALAMh«/ ) 

0 

c 

STOP 

0 

g/VD 
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SUbROLi f INE TRAwbKN (AAb»bS »bb»(NO*A»Pl »Hi IL/Eli) 

DIMENSION AXS(eO)» PS(bO)» HS(«0) 

IF (X .NE. XXS(MIN)) GO TC 75b 
Pi = PS(MIN) 

HI = hb(MIN) 

GO TO 756 

755 FACTOR = 1. 

MAX = MIN ♦ I DEG 
DO 701 J=MINtMAX. 

701 FACTOR = FACTO«*(X-XXS(J) ) 

c evaluate interpolating PCLYnOMIAL . 

PI = 0. 
hi = 0. 

DC 702 I=MIN»MAX 

TERMP = PS(I)*FACTOR/(X-XXS(I) ) 

TERMH = HS(I)<^>-ACTOR/(X-XXS(I) ) 

DO 703 J=MIN*MAX 

IF(I . NE. J) C30 TO 757 

GO TO 703 

757 TERMP = TERMP/ (XAS ( I ) -XXS ( J) ) 

TEkMH = TERMH/ (XXS ( 1 )-XXS (J) ) 

703 CONTINUE 

Pi = Pi ♦ TERMP 
HI = hi ♦ TERMH 

702 CONTINUE 

756 RETURN 
C 

EK'D 

SUdkOUllNE INTERPN ( Y »hlEF »Epi ) 

C 

COMMON BSTAR*X3,C1 
COMMON ALPHA»ALAMDA»GINV 

COMMON KX.NX1>NY.NYI,01»CX (36) 1 OY ( lA) »DEL (29. 15) » HTEp(29,15)» 
C hlHTEi- (29, 15) ♦ PK29), hl(29)» x(29) 

COMMON A(^9»1A»1h) » B ( 29 ♦ i *♦ » 1 5 ) ♦ C(29»14*1A). R(29»1A) 

COMMON Phl(29»iA)»PlPhl(2S»lA)»HT(29,ib) » NPR INT ♦ NPR INT2 
common PH10P(2‘:') 

DIMENSION Y (15) »H1EP (29) , EPl (29) 

C 

DO 406 I=1.NX 
DO 407 J=1»NY 

XYI = X(I)**2 ♦ Y(J)«*2 - 1. 

IF(XYI .GE. 0.) GO TO 451 
DEL(I.J) = C1*XYI 
6C TO -+52 

451 DEL(I.J) = 0. 

452 HT(I.J) = Hid) ♦ DELdfiw) 

HTEP(IO) = HT(I»J)**3*£F1 (I) 

H1HTEP(I*J) = HIEP(I) - HTEP(I.J) 

407 CONTINUE 
406 CONTINUE 
J=1 

PRINT 416 

PRINT 7» (I.OELd.J) »I = 1*NX) 
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Pt^INT ^17 

PPINT 7*(I« HTd»J) »I = 1»KX) 

IF(NPfilNT .EQ. 0) GO TC 490 
PPInT h18 

PRINT iO» ( (HTfc-PdtJ) »J=1*NY)»1 = 1»NX) 

PRINT 419 

PRINT *0» ( (HlriTEPd* J) t'v-=l»NY) » I = i»NX) 

C 

7 FORMAT (8 (X» I2»X»EI3,6) ) 

10 FORMAT (Xd5E9.i) 

416 FORMAT T>£L«/ ) 

417 FORMAT I //<&y, * H I «»/ ) 

418 format ( //6X* -^H I EP</) 

419 format (// 6X.*hiPTEP*/) 

C 

490 RETURN 

C 

END 

SUbROU-flNE ABCk 
C 

COMMON BSTAR»XJ»Ci 

common alpha*alamca»ginv 

common NX,NX1 ,i\Y*.\Yi »D1 .CX (28) » DY ( 1 4 ) , DEL ( 29 d5 ) » HTEP(29*15)» 

C PlHTEP(29d5) * r'l(.;9), l-l(29)» X(29) 
common A(29»14*1a)« B { 2 “7 » i ♦ 1 5 ) * C(29»14*14)» k(29»14) 

COMMON Phi (29»i*i) .PlPhl (EStlHJ »HT(29»lb) * NPR INT ♦ NPH INT2 
COMMON PhIOP (29) 

C 

DO 501 I=1»NX 
00 501 J=1»NY1 
DO 501 K=1»NY1 
A(I,J,K) = 0. 

B(I,J»K) = 0. 

C(1*J»K) = 0. 

501 CONTINUE 

DO 502 I=1»NX.NX1 
DO 502 J=itNYl 
P(I»J*u) = 1. 

R(I*J) = 0, 

502 CONTINUE 

DO 503 1=2«NX1 
J = 1 

Ad,J»J) = DY(J)/DX(I-i)« (hTEPd»J) +HTEP(I-1»J) ♦ALPHA*HTEP(1-1»J) 
A « (PI ( i )-Pl ( I-U ) ) 

ed,J.U) = -DY ( J) *( (riTEP ( I *i > J) -t-HTEP ( I » J) ) /DX d ) ♦ (HTEP ( I» J) ♦ 

B HTEP(i-l»J) )/OX(I-i ) + ALPhA<»hTEP(I»J)*( (PI (I*i)-Pl (I) )/Dxm- 
R (PI (I)-Pl d-1) )/DX(I-i) ) ) - GiNV*(OX(I-l) +DX(I) )*(HTEP(I*J^1)^ 

B PTEP ( 1 • J) ) /DY ( J) 

Bd.J»J»l) = GI NV*^(0X(1-1)+0X (!))*( HTEH(dJ+l)+hTEP(I»J) )/DY ( J) 
Cd«J»J) = DY(J)/D>. (I)*(PTEP(dl» J)*HTEP(i»J)-ALPhA*HTEP(I*l»J)* 

C (PI <i*l)-Pi (I) ) ) 

R(I. J) = uY( J) (hihTEPd + itj) +H1PTEP( I» J) ) * (PI ( 1 ♦ 1 ) -Pi ( I ) ) /OX (I ) - 

R (HlHTtP ( 1 , J) ♦rtlhTEr ( I-l ) *(P1 (I)-Pl ( I-i ) )/0X ( I-l ) +ALAMDA* 

R (DEL(i*l,J)-DtL(l-l,J))) 
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c 


DDX = (DX(I-1)+DX(1))/^J. 

DO 504 J=2»NY1 

DCY = (DY ( J-1 ) +DY ( J) )/<i. 

A < I » J» J) = DDY/DX <1-1 ) * (l-TEP ( I » J) +HTEP ( I-i ♦ J) •*-ALPHA*HTEP ( 1-1 f J)* 

A (PI ( I >-Pl( I-l ) ) ) 

b(I»J»J-l) = GINV*0 DX/DV (w-1)»(HTEP(I«J)+HTEP(1*J-1)> 

8(1. J*J) = -DDV*( (riT£t-(I*l*J)+rtTEP(I»J) )/DX(I) +(HTEP(I«J)* 

B HTEP(I-1 *J) )/L'X(I-l) + ALPHA<»MTEP ( !♦ J)*( (PI <1) )/DX(I)- 

F (PI (I)-Pl (I-l) )/ax(I-i) ) ) - GINV*DDX*( (HTEP(I»J*1)*HTEP(I»J))/ 

B-DY ( J) ♦ (HTEPd » J) (I»w-1) )/DY ( J-1) ) 

B(1.J*J+1) = GINV«DDX/(JY(w)*(mTeP(I»J + 1)+hTEP(I»J)) 

C(I»J»J) = DDY/DX(l)*(tiTEF (l*i*J)+HTEP(I»J) - ALP8A*HTEP(I*lfJ)* 

C (P1(I*1)-P1(D)) 

R(I.J) = DDY* ( (H1HTEP( I + l » J) ♦HlHTEP( I« J) ) « (Pi (I^i ^"F’l ( I) )/OX(I) - 

R (MlHTtP ( 1 » J) ♦^^lUTEP ( I-l rw ) ) * (Pi ( I ) -PI ( I-l ) )/0)< ' 1 ♦ ALAMOA* 

R (UEL(l+l,J)-OLL (I-l.J) ) ) 

50** continue 
503 CONTINUE. 

C 

IF (NPPINT .Eu, 0) GO TC SSiQ 
PRINT 821 
DC 505 I=1«NX 
PRINT 525 , I 

505 PRINT i 1 < ( ( A ( I » J * K ) ,K= 1 t 1 ) * J = 1 » Ny 1 ) 

PRINT 522 

DO 506 1=1 »NX 
PRINT b25*I 

506 PRINT 11 • ( (P { I « JiK > ,K = 1 ,N Y 1 ) , J=1 ,NY1 ) 

PRINT ,623 

DO 507 1=1. NX 
PRINT 625.1 

507 PRINT X 1 ♦ ( (C ( I * J »K ) .N= 1 »N 1 1 ) » J=1 «NY1 ) 

PRINT 624 

PRINT 11. ( (P(I»J) .j=l,NYl> ,i=l.NX) 

C 

11 format (X, l'+E9.«:) 

521 format (////6X.<'A(I,j»K)6/) 

522 FCRmaT(////6x.<»B(I.J»K)*/) 

523 FOkmaT (////GX.^'C ( I . J.K ) */) 

52.. format <////6X.*^R ( I *J) */> 

525 FORMAT (/6X.«I=*» 13) 

590 RETURN 

C 

Evp — 
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SOhROoHNl TEHr-hi 
C 

COMMON BSTAR»XJ»C1 
COMMON ALPHA*ALAM0A»GINV 

COMMON NX»NXl»tMY»NYl»DI»CX (28) ♦ DY ( 14) » DEL (29» 15) * HTEP(29fl5)» 
C hlHTEM29»15) ♦ Pl(29)» 1-1(29) ♦ X(29) 

COMMON' A(29*14»14)» B(29»14»15)» C(29*14*14)» R(29*14) 
common phi (29,1h) .PIPHI (29.14) .HT (29, 15) .NPRINT.NPRINT2 
C 

DIMENSION 1(29.14,14). E(30,14.W). F(30.l4). RAF(l4) 
dimension TT(1^.14), TI(1'<.14). IR(14), E£(14) 

C 

DO 601 I=1.NX.NX1 
00 602 J=1.NY1 
00 603 K=1.NY1 
T(I.J.K) = 0. 

E(1*1.J.K) = 0. 

603 CONTINUE 
T(I.J.J) = 1, 

F(I*1.J) = 0. 

602 CONTINUE 
601 CONTINUE 
C 

DO 604 1=2. NXl 
DO 605 J=1.NY1 
DO 605 K=1.NY1 
TT(J.K) = B(I.J.K) 

DO 607 L=1.NY1 

607 TT(J.K) = TT(J.K) ♦ A ( I .u »L) *E ( I .L .K) 

605 CONTINUE 

C 

C call library SUBROUTINE FCR MATRIX INVERSION 
CALL Mi (TT.NYl»NYl,DET.ICcT.TI»IR.IEP.EE) 

DO 608 J=l,NYl 
DO 608 K=1,NY1 

608 T (I.J.K) = TI (J.K) 

DO 609 J=1.NY1 
RAF(J) = R(I.J) 

DO 609 L=1»NY1 

609 RAE(J) = kAF(J) - A(I .J.L)*F (I’D 
DO 610 J=1.NY1 

Fd + l.J) = 0. 

DO 610 L=1.NY1 

610 F(I.l.J) = F(I + 1,J) -» T (1,J,L)*RAF(L) 

DO 611 J=l.NYl 

DC 611 K=1.NY1 
E(I*l»u.K) = 0. 

DO 612 L=1.NY1 

612 E(I*1«J«K) = E(I*1.J»k) “ T { I . J.L) ( I »L ,K ) 

611 CCNTlNoE 

604 CONTINUE 


144 



c 

NX2 = NX ♦1 

IF (NFHINT ,EU, 0) GO TC 6S»0 
PRINT t>61 
DO 631 I=1»NX 
PRINT 666*1 

631 PRINT i-l* (<T(I*J*K) ,K=l*NYl) *J=1,NY1) 

PRINT 062 

DO 632 I=2*NX2 
PRINT 066*1 

632 PRINT 11* ( (E(I*J*K) *K=l*Ntl) , j=1,NY1) 

PRINT 063 , . , 

DO 633 I=2*NX2 

633 PRINT il* (F(I*J) *J=1*NY1) 

C 

bRO CONTINUE 
I = NX 

DO 621 J=1*NY1 
Pbl (I *J) = F (I+l* J) 

621 PlPHi(i,j) = Pi(l) + PHI(I*J) 

DO 622 II=l*NXl 

I = NX-II 

DO 623 J=1*NY1 

Phi ( 1 * J) = F (I-i-l* J) 

DO 62^ K=1*NY1 

624 PHI(I*J) = PHld.J) + E(I+i*J*K>*PHI(I + l*K) 
623 PlPHl(i.J) = Pi(I) * PhI(I»J) 

622 CONTINUE 

IF (NPHINT2 ,EQ. 0) GO TC 691 
PRINT 064 
DO 625 I=1*NX 
PRINT 066*1 

625 PRINT /*(J*PH11I»J)*J=1*NY1) 

PRINT 065 

DO 626 I=1*NX 
PRINT 066*1 

626 PRINT 7* (J*PlPhI n»J) ,j=l,NYl) 

691 CONTINUE 

C 

7 FORMAT (8 (X* I2*X*E13.6) ) 

11 FORMAT (X*14E9. 2) 

661 F0RMAT(////6X**T(I*J*K)*/) 

662 format (////6X**E ( I *J*M*/) 

663 format l////6X*-F (I *J)«/) 

664 FORMAT (////6X**PHI ( 1*J)»/) 

665 FORMAT <////6X» =»PlPhl ( I * J)«/) 

666 F0RMATt6X*«I=*«lJ) 

C 

RETURN 

C 

END 


145 



OOOOOObOuOOOuOOOOOOOuO £nD Of ptCOKU 
112 110 

0. 03125 -0,25 -5. 

1 . 10 . 

0.3 

7.5^3b60E-13 0. 

' 0.003 100. 

0.000001 

61 69 29 29 15 2 

0.5 0.5 0.5 0.5 0.25 0,25 0.25 0.25 

0.25 0.c5 0.125 C.12'5 0.0625 0.0625 0.03125 0.03125 

0.015625 0.015625 0.00781250.0075125 U. 0078125 0.0078125 0.0078125 0.0078125 

0.0O78125 0.0078125 0.015625 C. 015625 0.03125 0,03125 0.0625 0.0625 

0,03125 0.03125 0.03125 C.03ic5 0.03125 0. 03125 0.015625 0,01562 

0.0078125 0.0078125 0.0078125 C.007dl25 0. 0078125 0.0078125 0.0078125 0.00781 

0,015625 0.015625 0.03125 C.031c5 0.03125 0.03125 0.03125 0.03125 

0.0625 0.0625 0.0625 0.0625 0,125 0.125 0,125 0.125 

0.125 0.125 0.125 0.125 

0, /21^50E-06 1.258592E-05 2.681879E-05 5.580325E-05 

8.104526E-05 1.220191E-0‘* 1.939579E-04 3.351647E-04 6.6'*4889E-04 

1.7068A0E-03 3.v5085‘»E-03 6.8679H3E-03 1.135745E-02 2.187715E-02 

3.203593E-02 4.':'93023E-02 6.344260E-02 b,l63219E-02 9.28**604E-02 

1.056967E-01 1.202896E-01 l.3b6692E-01 i.547954E-0i 1.745317E-01 

1 .956‘*41E-01 2.i7817tiE-01 2.o-;669yE-01 3.073501E-01 3.900927E-01 

4.590188E-01 5 . 079945 E-U 1 6.512e9‘+E-01 fa.86b959E-01 7.187712E-01 

7.480454E-01 7.748653E-01 7.995005E-01 8.221692E-01 3.328140E-01 

8.430270E-01 8.‘*79759£-01 8.5<;8227E-01 b.57569<4E-01 8.622175E-01 

8.667687E-01 8. /12246E-01 8.755866E-01 8.798561E-01 8.831239E-0i 

8.960368E-01 9.i08461E-01 9.243499E-01 9.366079E-01 9.476705E-01 

9.575759E-01 9.o63579£-Cl 9.8u76l8E-01 9.9O6762E-01 9.974024E-01 

r.000000E»00 

1.600004E+03 1.265680E*03 9.bd2892E*02 7.080937E+02 4.85‘+906E*02 

3,884724E*02 3.011240E+02 2,23608bE»02 1.561616E+02 9.914780E+01 

5.318671E+01 3.H69588E+01 l,954757t+0i 1.337916£+01 8.292149E+00 

b.207B85E + 00 4.‘+7262<r£ + 00 3.747008E + 00 3.122434E + 00 2.84dy27E*00 

2.601504E*00 2.j79973E*C0 2.1o382oE+ 00 2.012162E+OC 1.86370eE*00 

1.73b796£*00 1.O29467E+00 1 .464292E+00 1.350633E+00 1.219590E^00 

1.152006E*00 l.u95440E*uO 1.07l39bE*00 i.062275E*00 1.055176E+00 

1.049277£*00 1,u44232E*00 1,0j 9955E*00 1 .0361 aiE+00 1.034440E+00 

1.032793E+00 1.03d001E+00 1.0Jl22«£+00 1.030474£*00 1.U29738E+00 

1.029019E+00 1.u2b316E*00 1.J27628E*00 1.026954E*00 1.O25648E+00 

1.024392E*00 1 . o22o3'<E + 00 1 . u i 9b‘+<iE +u0 1.017787E + 00 1.015650E*00 

1.014002E+00 l.ul2247£*00 1.0U9023E*00 1.0C5711E+00 1.003650E*00 
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